{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Guided Inquiry (Rowe & Neshyba, 2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimal inverse retrieval of cloud properties "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**OBJECTIVE**: Build intuition and skill in retrieving atmospheric properties from remotely-sensed data\n",
    "\n",
    "**SELF-ASSESSMENTS**:\n",
    "\n",
    "1.\tDefine the terms _Planck blackbody function_ and _greybody model_.\n",
    "2.  Sketch a greybody spectrum, and explain how it would change when temperature and emissivity are increased or decreased.\n",
    "3.  Define and explain relationships between _observed quantities_, _quantities of interest_, and _forward models_ in the context of inverse retrieval, with examples.\n",
    "4.  Describe complicating factors that inverse retrieval methods typically encounter.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "It is often the case in atmospheric science that one would like to know some property of the atmosphere that is difficult to measure directly. For example, perhaps you'd like to know the temperature at a certain altitude, but you don't have an airplane to carry a thermometer to the necessary height. In such situations, one can fall back on indirect methods. One indirect method is called _inversion_: you remotely measure a quantity that _depends_ on the property you want, and from that you try to retrieve the property you are really interested in. Inverse retrieval is the focus of this exercise.  \n",
    "\n",
    "At the heart of any inversion problem is something called the _forward model_, designated _f_. If we use _x_ to symbolize the property or properties we're interested in (we'll call these the _quantities of interest_), and _y_ to symbolize the remotely-sensed observations (call these the _observed quantities_), then the forward model is a function that lets us calculate _y_ from _x_. Mathematically, we'd write this as _y=f(x)_. The forward model is generally an approximation to real atmospheric processes, but we hope that it is a good enough replica of what goes on in the real world that we can get useful information using it.  \n",
    "\n",
    "In fact, what we need is actually the inverse of the forward model: we need a way to calculate quantities of interest from observed quantities. Mathematically, we'd say we want $x=f^{-1}(y)$. Inversion is not always easy because of several complicating factors. First, there might be multiple quantities of interest (_x_) that are consistent with the observed quantities (_y_); then you have a non-unique inversion problem on your hands. Then there is the possibility that since the forward model is not an exact replica of what goes on in the real world, it could have systematic biases built in. And what if your observations are noisy? While there are multiple inverse retrieval algorithms available for accomplishing this, the one you'll learn about here is a powerful one that we'll call the _optimal inverse retrieval_ (see CD Rodgers, Inverse Methods for Atmospheric Sounding, 2000 (World Scientific Publishing)). This algorithm is optimal in the sense that it tries to fit the signal in the observed quantities, but not the biases and noise.\n",
    "\n",
    "Here, we'll focus on inverse retrieval of cloud properties as a concrete example. Suppose the cloud properties we are interested in are the cloud's temperature and thickness, but we can't measure those quantities directly. Instead, we have observations of the infrared radiance coming from the cloud at a set of frequencies. So the quantities of interest (_x_) are the cloud temperature and its thickness, and the observed quantities (_y_) are the cloud's infrared radiance spectrum, as measured by an instrument on the ground looking up at the cloud. We'll need a forward model that relates these, of course. A simple one is called the _greybody model_; in the greybody model, the cloud is modeled as a Planck blackbody function multiplied by a factor that accounts for the cloud thickness. Mathematically, we write the greybody model as\n",
    "\n",
    "<p style='text-align: center;'>\n",
    "$ y = \\epsilon B(\\nu,T) \\qquad $(1)\n",
    "</p>\n",
    "\n",
    "where $B(\\nu,T)$ is the Planck blackbody function, $\\nu$ is the frequency of infrared light, _T_ is the cloud's temperature, and $\\epsilon$ is called the _emissivity_ of the cloud. The emissivity is a measure of how well the cloud emits, which in turn depends on the cloud's thickness: $\\epsilon\\approx0$ if the cloud is extremely thin (or if there is no cloud at all), and $\\epsilon\\approx1$ for an extremely thick cloud. In this equation, we have as many values of $B$ (and therefore of _y_) as we have frequencies $\\nu$. \n",
    "\n",
    "The forward model used here involves a variety of approximations. For example, in real clouds $\\epsilon$ changes with frequency, and _y_ is influenced by the atmosphere between the ground and the cloud. But we will ignore those complications here.\n",
    "\n",
    "The inverse retrieval solution to this problem is set up using matrices and vectors: in this case, _x_ is a vector with the values _T_ and $\\epsilon$, and _y_ is a vector with the values of observed radiance at a set of frequencies. The algorithm is  iterative, which means you might need to apply the algorithm a few times before converging on a good answer. Because it's iterative, we have to talk about _current solutions_ and _next solutions_. We'll designate the current solution as $x_n$, and the next solution as $x_{n+1}$. Then we have\n",
    "\n",
    "<p style='text-align: center;'>\n",
    "$ x_{n+1} = x_a + (S_a^{-1} + K_n^tS_\\epsilon^{-1}K_n)^{-1}K_n^tS_\\epsilon^{-1}[y-f(x_n)+K_n(x_n-x_a)]  \\qquad $(2)\n",
    "</p>\n",
    "\n",
    "where $x_a$ is the a priori state of the desired quantities (state vector), $S_a$ is the variance (error) of $x_a$, $K_n$ is weighting function for the particular inversion problem, $S_\\epsilon$ is the variance (error) in the measurement, y is the measurement, and f designates the forward model [where y = f(x)].\n",
    "\n",
    "This equation is the same as eq 5.9 in Rodgers (2001), Chapter 5 on <i>Optimal Nonlinear Inverse Methods</i>.\n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Created on Thu Apr  2 09:05:09 2015\n",
    "authors: Penny Rowe and Steven Neshyba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Resources needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "solve = np.linalg.solve \n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Define the blackbody function as a function of wavenumber"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plancknu(nu_icm_in,T):\n",
    "    import copy\n",
    "    nu_icm = copy.deepcopy(nu_icm_in)\n",
    "\n",
    "    # spectral Planck function as function of wavenumbers (cm^-1)\n",
    "\n",
    "    # [h]    = J*s\n",
    "    # [c]    = m/s\n",
    "    # [cbar] = cm/s\n",
    "    # [k]    = J*K-1\n",
    "    # [B]    = cm*s-1*J*s*cm3*s-3*cm-3*m-2*s2\n",
    "    # [B]    = W m-2 cm\n",
    "    #\n",
    "    h    = 6.62606896e-34             # J s;  CODATA 2006\n",
    "    c    = 2.99792458e8               # m/s;  NIST\n",
    "    k    = 1.3806504e-23              # J K-1; CODATA 2006\n",
    "    cbar = 100*c                      # cm/s\n",
    "\n",
    "    indzero = np.where(nu_icm==0)            # avoid divide-by-zero  \n",
    "    nu_icm[indzero]=.1\n",
    "  \n",
    "    top        = 2 * h * cbar**3 * nu_icm**3\n",
    "    bottom     = c**2*  ( np.exp(h*cbar*nu_icm/(k*T))-1 )\n",
    "    f          = cbar*top/bottom \n",
    "    f[indzero] = 0\n",
    "    return f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Define the forward model, with units $mW \\space / \\space (m^2 \\space sr \\space cm^{-1})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def greybody(nu,X):\n",
    "    T   = X[0]\n",
    "    Eps = X[1]\n",
    "    y = 1e3*plancknu(nu,T)  * Eps\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Define a set of wavenumbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nobs = 20\n",
    "nu = np.linspace(200,1500,Nobs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font face=\"georgia\">Pause for Analysis #1. Use graphics to get a sense of the forward model, by graphing the greybody radiance as a function of wavenumber for these values of temperature and emissivity, and then perturbing the temperature and emissivity a little. For example, how does the maximum radiance change when the temperature is increased by 100 degrees? When the emissivity is increased to 1? How does the _frequency_ at which the maximum radiance occurs change when these changes in temperature and emissivity are made? Make appropriate sketches to record your results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1235ebf98>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3yO1//H8dfJlp2QkCCJERp7xN6ColS3lrZGVLWlRcdP9eurvmiLUlpKbVqKqtJqaxSxV6gYsSMRITJljzvJ+f1xpa22Rsg9k/N8PPJIcst9nY8+9O1yrs85R0gpURRFUSyPlakLUBRFUR6OCnBFURQLpQJcURTFQqkAVxRFsVAqwBVFUSyUjTEHq1SpkgwICDDmkIqiKBbv2LFjSVJKr3++btQADwgIIDw83JhDKoqiWDwhRMydXldTKIqiKBZKBbiiKIqFUgGuKIpioVSAK4qiWCgV4IqiKBZKBbiiKIqFUgGuKIpioYzaB648nBMJJziecJyqzlUJcA2gukt1HG0dTV2WoigmpgLcjJ1IOMH8iPkcuH7gX79W2bEy/q7+f34EuAbg7+pPVZeq2FrZmqBaRVGMTQW4GTqZeJIvT3zJ/uv78bD3YGzzsfSr3Y/E7ESi06OJSY8hJj2G6PRotsVsIy0v7c/3WgtrqrlU+yvcXfzxd/OnjkcdPB08Tfi7UhRF31SAm5FTiaf4MuJL9sXtw8PegzHNx/B83ef/nC7xdPCkrmfdf73vVu4tYjKKQz3tr4A/Gn+UnIIcAGytbBneaDihDUKxtVZ36IpSFqgANwOnk07z5Ykv2Ru3F3d7d0Y3G80Lj7xQ4nludwd33B3caezV+G+vSylJyE4gOj2a9RfWM+/EPLbFbON/bf9Hg0oNDPFbURTFiIQxz8QMDg6WajOrv5xJOsOXEV+y59oe3OzdGFx/MC888gJOtk4GGW/X1V1MOTyFpJwkXgx6kTeavKEehiqKBRBCHJNSBv/zdXUHbgJnks8w/8R8dl/bjZu9G282fZMBQQMMFtx/6OLXheAqwXx27DNWRq5kx9UdfNj2Q1r7tDbouIqiGIa6AzeiyORI5p+YT9i1MFztXBlUfxADHhmAs52z0Ws5Gn+USQcnEZMew5O1n+Tt4Ldxs3czeh2Kotzf3e7AVYAbwbWMa0w/Op1dsbtwsXNhUL1BDAwaaJLgvl1uQS4LIhaw/MxyPBw8GN9qPN39u5u0JkVR/k0FuIlk5Gcw4OcBJOUkMai+Ftwudi6mLutvziafZeKBiZxNOUuIXwgftPoAL8d/Hf6hKIqJ3C3A1VJ6AyqSRYzfN57YjFg+7/o5IxqPMLvwBgiqGMSqx1Yxutlo9sXto9/Gfmy4uAFj/uWuKMqDUwFuQAtPLiQsNox3gt+hRZUWpi7nnmytbAltGMr6vuup41mHiQcm8sq2V4hNjzV1aYqi3IUKcAPZc20PX574ksdqPsbAoIGmLqfEAtwCWProUia0nsCZ5DM89eNTrDizgoKiAlOXpijKP6gAN4Cr6VcZt2ccdT3rMrHNRIQQpi7pgVgJK56r+xwb+22ktW9rPg3/lEG/DiIlN8XUpSmKchsV4HqWrcvmrV1vYWVlxWedP6OCTQVTl/TQKjtV5vMunzO943TOp54ndGsoyTnJpi5LUZRiKsD1SErJxAMTiUqLYnqH6VRzqWbqkkpNCEGvGr2YFzKPuMw4QreGkpSTZOqyFEVBBbherYxcyZboLYxqOoq2Vduauhy9auXTinkh87iedZ0hW4aQkJ1g6pIUpdxTAa4nh28cZtaxWXT3705og1BTl2MQLaq0YEG3BSRkJzB061Dis+JNXZKilGslCnAhxBghxBkhxGkhxLdCCAchhKcQYrsQ4mLxZw9DF2uubmTe4N3d7xLgGsDkdpMt7qHlg2hWuRlfdf+KpJwkhmwZwo3MG6YuSVHKrfsGuBCiKvAmECylbABYA88D44AdUspAYEfx9+VObkEuo8NGoyvSMbvLbINvSGUOmng3YWH3haTlpTFk6xDiMuNMXZKilEslnUKxASoIIWwAR+A60A9YUfzrK4An9F+eeZNSMuXQFCKTI/mo/UfUcKth6pKMppFXIxb1WERGfgZDtgwhNkMt+FEUY7tvgEsp44BPgavADSBNSrkNqCylvFH8MzcA7zu9XwgxXAgRLoQIT0xM1F/lZmDd+XVsuryJVxu9She/LqYux+jqV6rP4h6LyS7IZsiWIVxNv2rqkhSlXCnJFIoH2t12DcAXcBJCvFjSAaSUC6WUwVLKYC+vsrNB0omEE3xy9BM6VO3A601eN3U5JhNUMYglPZaQX5jPkC1DiE6LNnVJilJulGQKpRtwRUqZKKXUARuAtsBNIYQPQPHnctNXlpidyJiwMfg4+fBxh4+xEuW7maeuZ12WPLqEAlnAkK1DiEqLMnVJilIulCR5rgKthRCOQmuvCAHOAj8Cg4p/ZhCwyTAlmhddoY6xYWPJ0mUxu8tsdQhCsUCPQJY+uhQpJUO2DOFS6iVTl6QoZV5J5sAPA+uB48Cp4vcsBD4BugshLgLdi78v86YdncaJxBNMajuJOh51TF2OWanlXotlPZdhLawZunUo51POm7okRSnTSvRvfynlRCnlI1LKBlLKl6SUeVLKZClliJQysPhzmd/paOOljaw9v5ZB9QbRq0YvU5djlmq41WBZz2XYWtsybNswzqWcM3VJilJmle/J2wdwJukMkw9OplWVVoxuPtrU5Zg1f1d/lj+6HAcbB0K3hnIm+YypS1KUMkkdqVYCKbkp9N/cH4FgTZ81eDp4mrokw5ASdNmQk1r8ceuvr4WAgPbgWbPEl7uWcY3QraFk5GfwVfevaOjV0IDFK0rZdbcj1WxMUYylmRk+k5ScFFb2WmmZ4Z2ZCFFhtwVzKuTeunNQF+nufS2PGlA7BGp1hYAO4OB61x+t5lKNZT2XEbo1lOHbh7Oy10oCPQL1+3tTlHJM3YHfR2x6LH039mVg0EDebfGuqct5MOnXYf/ncGw5FOT89bqdC1TwgAruxR8ef3043P79bV/nZ2t/CVzeAVf2gi4LrGygWkuo3RVqhYBPE7D696zcjcwbDPhlAHZWdqx6bBWVKlQy2n8CRSkL1Kn0D2nigYlsvryZLU9vsZyT2lOjYd9sOLEKigqh8fPQ8hVwqw4ObmBtW7rrF+RD7GEtzC/tgPiT2usVPKFWFy3Ma3UFV58/33I66TRDtgyhjmcdlvRYgoONQ+lqUJRyRAX4Q4jLjKPPhj48V/c53m/1vqnLub+kS7B3JpxcC1bW0PRFaDcaPPwNO25mIkTtgss7tY/Mm9rr3vW0IK8dAn5t+e36PsaEjaFXQC+mdZxWpndtVBR9UnPgD2HpqaUIIRjSYIipS7m3m2e04D7zA1jbQ6tXoe0ocPU1zvjOXtDoOe1DSrh5WgvySzvgyEI4OBfsXenW7i3eavw6cyK+xN/NnzeavGGc+hSljFIBfhfxWfFsuLSBp2o/RRWnKqYu586u/w57PoVzm8HOGdq+CW3eAOc77itmHEJAlYbaR7u3ID8LYg5A+FLYOZlQ5yrEBDZlQcQC/Fz86Furr+lqVRQLpwL8LpaeXgoSQhua4ek6Vw/Dnhlwabs2p91pnHbX7WiGHTJ2ThDYXfuIOYjY/l/++/uvxFUPYOL+CVR18qVZleamrlJRLJJayHMHidmJfH/he/rV7oevs5GmIe5HSriyB5b3gaU94PpxCJkIo09Dl/fNM7z/yb8NhG7Dtv8qPsu2wTc/l9FbhhJ7/kdTV6YoFkkF+B0sO7OMQlloPnffN07C0kdhRV9IugiPfgSjT0GHsffswzZLQkBQH9xeP8y8eq9SJAt5Y897pK9+DhLUsntFeRAqwP8hOSeZ785/x2M1H6O6S3VTlwMxB2H5Y5AaA4/NhLcitHluOws/us3aBv+2Y/is+3xi7ewZm3kS3fw2sGmk1r+uKMp9qQD/hxWRK8gvyueVhq+YuhStk+PrJ8G5MryyE1oMA9uy1T/doloHJrabxGF7W6YGtUVGrIHPm8Fvk7QVooqi3JUK8Nuk5qay5twaetXoRYBbgGmLObsZVveHirVhyK/gVtW09RjQE7WfYFjDYXyfc5WVvT6AoL6wbxZ83gQOzoOCPFOXqChmSQX4bb6O/JrcglyGNxxu2kJOroN1L4NPYxj8k9ZnXcaNajqK7v7dmRm5jJ0tXoBX94BvU9g6Hr4IhpPfaQ9yFUX5kwrwYml5aaw+t5oeAT2o6V7yHff0LnwpbBgO/m3hpY3aPiTlgJWwYmr7qdSvWJ9xe8cRaWcLL/1Q/N/AHTYM0/5Sy0k1damKYjZUgBdbdXYVWboshjcy4d33/s9h8xgI7AEDvwN7Z9PVYgIVbCrwRcgXuNm7MWrHKG5m3dT2Vhm+G7r/D87/AvPbaw92FUVRAQ6QkZ/BN5HfEOIXYppj0qSEXR/B9glQ/0no/w3YVjB+HWagUoVKzO06l0xdJqN2jiJbl63tcNjuLQjdpm3Etbw3hE3TNupSlHJMBTiw+uxqMnQZvNroVeMPLiVs/QB2T9M2n3p6CdjYGb8OM1LXsy4zOs3gfOp5xu0dR+EfQV21uTY33uAZCPtI64tPizNtsYpiQuU+wLN0WXx99ms6V+tMUMUg4w5eVAg/vQWH5kGrEdD3C20XQYWO1TryXov32BW7i9nHZ//1Cw6u8PQiePIruH4CFrTTOnYUpRwq9wG+5twa0vLSeLWxke++C3Xww6twfAV0eAd6fnLHwxDKswGPDKB/3f4sP7Oc7y98//dfbPw8jNgL7n6wdiD8/Dbocu58IUUpo8p1YmTrsllxZgXtqrajQaUGxhtYlwvrBsGp77T9TEImaEvMlb8RQjCu5TjaVW3HlENTCI//x17yFWtB6G/QZiQcXQyLukLCWdMUqygmUK4D/LsL35Gal8qIRiOMN2h+FnzbH87/DL0/1fYzUe7KxsqGGR1nUM2lGm/vfpv4rPh//IAdPDoVBn4PWYmwsAuEL1M940q5UG4DPLcgl2Wnl9HapzVNvJsYZ9CcW9rS+Ct74In52jFnyn252Lkwu8tscgtyGRs2lrzCO6zMDOwGI/aDX2vYPFr1jCvlQrkN8O8vfk9ybrLxOk+ykrSuibjj8MwyaDLAOOOWEbXcazG1/VROJZ1i6qGp3PEoQJfK8OIG1TOulBvlMsDzCvNYemopwZWDCa7yr2Pm9C8jHpb1hqQL8MK3UP8Jw49ZBnXz78YrDV/hh0s/8N2F7+78Q6pnXClHymWAb7y4kYScBEY0NsLcd2GB9sAy7RoMXK+dTKM8tDeavEH7qu35+MjH/J7w+91/8E4945mJxitUUYyg3AW4rlDH4tOLaerdlJZVWhp+wH2fQewh6PMZ1Ohg+PHKOGsra6Z1nIavky9jw8aSkJ1w9x/+o2f8iQXa1NXirnAz0njFKoqBlbsA33R5E/FZ8bza6FWEoVv3rh2DsI+hwdPaie2KXrjauTK7y2yydFmMDRtLfmH+vd/Q5AUY8gsU5MOSHnBxu3EKVRQDK1cBrivSsfjUYhpWakhb37aGHSwvEza8Ai4+8Ngs1eetZ4EegUxpN4WIxAg+OfLJ/d9QtZl2KIZnAKx+Dg4tUK2GisUrVwH+c9TPxGXGMaLxCMPffW8dDylR8NRX2naoit71COhBaINQvrvwHesvrL//G9yqwpAtUKcXbPk/+HmstiJWUSxUuQnwgqICFp1cRJBnEB2qGngu+uxmbYl8u7cgoL1hxyrnRjUdRVvftnx0+CMiEiPu/wZ7Z223x3ZvaXuvr3pGHd2mWKxyE+BbordwNeMqrzY28Nx3Rjz8OAqqNIIuHxhuHAXQHmpO7zgdb0dvxu4aS1JO0v3fZGWl9Yr3mwfR+2FJd0i+bPhiFUXPykWAFxYVsvDkQup41KFL9S6GG0hK2Pg66LLh6cXlfltYY3Gzd2NOlzlk6DJ4O+xtdCWdFmn6Iry8UVuCvzhEC3NFsSDlIsC3x2znStoVhjcajpUw4G/5yEK4vAN6TAGvuoYbR/mXup51mdR2EscTjjP96PSSvzGgPQzbAY6VYGU/+H2V4YpUFD0rFwG+5vwa/F396e5vwEU0CWdh2wTtOLQWwww3jnJXvWr0YnD9waw5v4YfLv5Q8jdWrAXDtkNAO9j0OmyfCEVFhitUUfSkzAd4bHosx24e44naTxju7rsgD74fBvYu2ryqahk0mbeavUUrn1ZMOTSF00mnS/7GCh7aStngobB/Nqx7Sds5UlHMWJkP8B+jfkQg6FOzj+EG2TkZbp7WwtvZ23DjKPf1x/azlSpUYvSu0STnJJf8zda2Ws9+z2naZlhLe6oj2xSzVqIAF0K4CyHWCyHOCSHOCiHaCCE8hRDbhRAXiz97GLrYB1Uki/jx0o+08W1DFacqhhkkajccmKvdudXtaZgxlAfi4eDB7C6zuZV3i3d2v4Ou6AF6vYWA1iPghbWQckU7JCLuuOGKVZRSKOkd+Bxgi5TyEaAxcBYYB+yQUgYCO4q/Nyvh8eFcz7pOv1r9DDNAdgr8MAIq1oYeUw0zhvJQgioGMbHNRMJvhjMrfNaDX6BODwjdCtZ22k6SZzbqv0hFKaX7BrgQwhXoCCwBkFLmSylvAf2AFcU/tgIwuz1SN13ehLOtM139uur/4lLC5jGQlaBtmGTnqP8xlFLpW6svLwa9yDdnv+Gnyz89+AUq19eW31dpCN8NgoPz9F+kopRCSe7AawKJwDIhxO9CiMVCCCegspTyBkDx5ztO/gohhgshwoUQ4YmJxtvOM0uXxfaY7fSs0RMHGwf9DxDxLURuhC7jwbep/q+v6MXY4LG0qNKCSQcncSb5zINfwNkLBv0EQY9r2yNs/UB1qChmoyQBbgM0A+ZLKZsCWTzAdImUcqGUMlhKGezl5fWQZT64bdHbyCnIMcz0ScoV+OVd8G8H7Ubr//qK3tha2fJpp0/xdPB88Ieaf17EAZ5dDi1egYNz4Yfh2s6GimJiJQnwa8A1KeXh4u/XowX6TSGED0Dx53tszGx8my5vIsA1gMZejfV74cIC+OFVENbw5AKwstbv9RW983TwZHaX2aTmpvL27rcf7KHmH6ysofcMCPkvnPoOVj8LeRn6L1ZRHsB9A1xKGQ/ECiH+WFoYAkQCPwKDil8bBGwySIUP4Y/e7361++l/35N9syD2MDw2E9z99HttxWDqVazHh20/5NjNY8w4OuPhLiIEdHgb+n0JV/bC8scg46Z+C1WUB2BTwp8bBawSQtgBUcAQtPBfJ4QIBa4CzxqmxAdnsN7va+EQ9gk0fBYamc1vVymhPjX7cC75HCsiVxDkGcSTgU8+3IWaDgQnL+3B5pLu8NIP2mpORTGyErURSilPFM9jN5JSPiGlTJVSJkspQ6SUgcWfUwxdbEkYrPf7jwMaXH2h96f6u65iVKObj6a1T2smH5rMycSTD3+hOj20h5t5GVqIxx3TX5GKUkJlbiWmwXq/t76vPbx8coE6oMGC/bFS09vRmzG7xpCYXYrOqGrBELod7JxgeR+4+Jv+ClWUEihzAW6Q3u9zv8DxldB+tDqgoQxwd3D/c/vZMWFj7n+m5r1Uqq2FeMVa8G1/OLFaf4Uqyn2UqQA3SO93QR5sGQfe9aDzeP1cUzG5up51+V+7/xGRGMHHRz4u3cVcqsDgX7S20o2vwd5Z6rxNxSjKVIAbpPf76BK4FQM9JqsDGsqYngE9CW0QyvoL61l3fl3pLubgqu1m2OAZ2DEJfv0/KCrUT6GKchcl7UKxCHrv/c65BXumQ80uULubfq6pmJVRTUdxLvUcHx/5mECPQJp6l2JVrY0dPLUInCvDoXmQeROe/EpbCKQoBlBm7sAN0vu9b5YW4t3/p5/rKWbH2sqaaR2m4evky5hdY7iZVcq+bisr6PmRdipT5Eb45ml1aLJiMGUmwPXe+30rFg4tgMbPg08j/VxTMUt/nKmZU5DDmLAx5BXmlf6ibUdpd+Oxh7TdDNOvl/6aivIPZSLADdL7vXOK9lmdLF8u1PaozUftP+JU0ikmH5yM1MdDyEbPwYB12jOUJT0g8ULpr6kotykTAa733u8bEXByLbR+Ddyr6+eaitkL8Q9hROMRbLq8iW/Pfaufi9YOgcGboSAXlvbQVvMqip6UiQDXa++3lLD9v9oZie3HlP56ikV5rfFrdK7WmelHp3M0/qh+LurbFIZuBQc3WNEXLm7Xz3WVcs/iA1zvvd+Xd0BUGHR6T624LIeshBUfd/gYP1c/3g57mxuZN/Rz4Yq1YOg27fPq/nBCT3f4Srlm8QGu197vokLY9l/wCIDg0NJfT7FIznbOzOkyB12Rjrd2vUVOQY5+LuxSWVvwE9AONo6A/XPUgh+lVCw+wPXa+x2xBhLOQMhEtWinnKvhVoNpHadxLuUckw5O0s9DTfhrwU/9J7Wpum3/USf8KA/NogNcr73f+dla50nV5tr/XEq517FaR0Y2HcnPUT+zMnKl/i5sYw9PL4WWw4tP+HlVnfCjPBSLXomp197vw/Mh4zo8vVjbuF9RgFcavsK5lHPMOjaLGm416Fito34ubGUFvaZrqzZ3TobsJHjua7B31s/1lXLBYu/A9dr7nZUEez+Dur21+UlFKSaEYEq7KdT1qMu7u9/lfMp5fV4cOr4DfT/XHpyv6Kv9WVSUErLYANdr7/fu6aDLhm6TSn8tpcxxtHXki65f4GzrzMidI0u3h/idNB8E/VdBQqS24Cc1Rr/XV8osiw1wvfV+J1+G8CXQ7GXwqqOf4pQyp7JTZeaGzCUtL41RO0eRrcvW7wCP9IaXN2lTKUt6QPxp/V5fKZMsMsD12vu9YxJY20Pn9/VTnFJmBVUMYnrH6UQmRzJ+33iKpJ67R/xaawt+hBUs6wXR+/R7faXMscgA11vvd+wRiNwE7d7UenQV5T46V+/Muy3eZcfVHcw+Nlv/A3gHQeg27ZCIr5+CyB/1P4ZSZlhkgOul91tK2DZB6wJoM1J/xSll3otBL9K/bn+WnVnG+gvr9T+Ae3XtTtynEax7WTtURFHuwOICXG+93+c2a1t9dn5ftW4pD0QIwbiW42hXtR1TDk3h4PWD+h/E0VObEw/sAT+PhV0fq1Wbyr9YXIDrpfe7UAe/fQiV6kLTl/RWm1J+2FjZ8GnHT6nhVoO3w97m8q3L+h/EzgmeXwVNBsLuT+CnN7U/u4pSzKICXG+938eWQ/Il6D4JrC16LZNiQs52zswLmYedtR1v7HiD5Jxk/Q9ibQv95kHHd+H4Svj2ecjL0P84ikWyqADXS+93XgaEfQL+7aFOT/0Vp5RLvs6+zA2ZS3JOMm/uepPcglz9DyIEdP0P9J0Dl3dpJ/xkxOt/HMXiWFSA66X3e/8crde2x//UknlFLxpUasBHHT7iZOJJJuyfoP/2wj80HwwvrNHWLizuBgnnDDOOYjEsJsD10vudfh0OzIUGT2ubVimKnnT3786Y5mPYEr2FeSfmGW6gOj1gyC9QmK8t+Lmy13BjKWbPYgJcL73fuz6CogLoOkF/hSlKsSH1h/BU4FMsPLmQTZc2GW4g3yYQul3rFf/mKThlgFZGxSJYTICXuvf7ZiScWKVt4elZQ7/FKQpae+F/Wv+HVlVa8eHBD/V3JNudePhD6Fao1gK+D4V9n6k2w3LIIgJcL73fv00EOxdt9zcLcjM9l00n4vh8x0VWHIhm04k49lxI5OS1W8SmZJORq9PfYQNKqdla2TKz80yqu1Rn9K7RRKdFG26wCh7w0g/alOBvH8LPb0NhgeHGU8yORfTQlbr3O2o3XNym7Tbo6Knf4vQsIT2XQ1dSOHg5mcNRyUQlZd33PbbWArcKdng42uLhaIf7H5+dtM8ejrb4ulegZQ1P7G2sjfC7KN/c7N2YFzKPgT8P5I0db7Cq9yrcHQx0vqqNPTy1GNyqaQ/o0+PgmaVaD7lS5glj3r0FBwfL8PDwB37f1uitnEo8xTstHuLuWUpY2hPSYmHUcbDVw8HHepSQkcvhqBQORSVzMCqZqEQtsF3sbWhZw5PWNSvSumZF6lRxJiO3gFvZOm5l55OarSM1O//Pr29l55Oa9cdrf33OL/yrI8LF3oauQd70alCFjnW8cLSziL+/LdaJhBOEbg2lQaUGLOqxCDtrAx/Td2QR/Poe+DSBAWvB2duw4ylGI4Q4JqUM/tfrlhDgpRK9D5Y/Br1mQKvhxh37DhIz8jh8JVkL7MvJXC4ObOc/A1sL7fq+blhbla7NUUpJjq6Q1GwdF+Iz2HI6nm2R8aRm63CwtaJzHW96NqhC1yBvXB1s9fHbU/7h1yu/8t6e9+hbsy9T208t/dF/93PuF1g/VAvvF7+HSoGGHU8xivIb4Cv7aQ8wR58E2wrGHbtY5PV0vj1ylUNRyVxMyATAyc6aFjU8aVN8h13f1xUba8M/kigoLOJIdApbTsez5XQ8CRl52FoL2tWuRK8GVeherwqeTupAZ336KuIr5p6Yy+D6gxnbfKzhQ/zaMVj9HMhCrW/cr7Vhx1MMrnwG+LVwWBwC3f8H7d4y3rjFsvMLmPPbRRbvu4K9jRUtArS76za1KtLASIF9L0VFkt9jb7Hl9A1+PR3PtdQcrAS0qlGRXg2r8Gj9KlR2Na8pJ0skpeTjIx/z7blvGdlkJK82ftXwg6ZEwTfPQNo1eGoh1H/C8GMqBlM+A3z183D1IIw5DfYuxhsX2H0hkQ9+OMW11Byeb1Gdcb0ewd3RfO9spZScuZ6u3ZmfiedS8b8Umvm506uBDz0bVKG6p6OJq7RcRbKICfsn8OPlHxnXchwDgwYaftCsZFjzgrbvfY8p0OYNtfrYQpW/AI8/DQvaQefx0Pn/jDMmkJSZx+TNkWw6cZ2aXk58/GRDWtWsaLTx9eVSgjZn/uvpeM5cT8dKQN/GvozsUpvAysb9y7CsKCgq4O2wt9kZu5Mp7abQr7YeznO9H10ObBgOZ3/UluL3mgE25nsjodxZqQNcCGENhANxUso+QghPYC0QAEQDz0kpU+91DaMG+HdD4OJ2GHNK65c1MCkl68Jj+eiXc+TkF/J6l1q81rlWmWjbu5qczaojMXx9MIYcXW9GfRMAACAASURBVCG9G/owqmttHqniaurSLE5+YT5v7HiDI/FHmNlpJt38uxl+0KIi2DkZ9s0C/3bw3NfgZHk3FeWZPgJ8LBAMuBYH+HQgRUr5iRBiHOAhpbznra7RAjzpIsxtAe1HQ7cPDT7c5cRMxm84xeErKbQM8OSjpxpS27vsHRKRkpXPkn1RrDgQQ2ZeAT3rV2FUSG3q+7qZujSLkq3LZvj24ZxJPsO8rvNoW7WtcQY+uQ42jdSOD3xhDVSub5xxlVIrVYALIaoBK4CpwNjiAD8PdJZS3hBC+ABhUsq697qO0QJ84+twegOMPgXOXgYbJq+gkAVhUczbdQkHWyvG9w7iueDqWJWy/c/c3crOZ+n+aJbtv0JGbgHdgirzZkhtGlUz0GKVMig9P52hW4YSkx7Dwh4Laerd1DgDXzsGawZAfiY8tQge6W2ccZVSKW2Arwc+BlyAd4oD/JaU0v22n0mVUv5rrkIIMRwYDuDn59c8JiamFL+NEkiNgc+bQstXoNc0gw1zNDqF9zec4lJCJn0b+zKhTxDeLuWrYyMtR8fKA9Es3neFtBwdnet68WZIIM38DD9lVRYk5SQxZMsQknOSWfLoEoIqBhln4PTrWohfPwEhE6D9WPVw08w9dIALIfoAvaWUrwshOvOAAX47o9yBbx6rnVzyVgS4VdX75dNydHzy6zm+PXKVqu4VmPJkA7rULd8r3jJydXx9KIZFe6JIzdbRIbASb4YE0iLAvLctMAc3Mm/w8paXyS/MZ1nPZdR0q2mcgXU5sOkNOP09NHwWHv/CZOsklPsrTYB/DLwEFAAOgCuwAWiBuU2hpN+AOY2h8fPw+Od6vbSUkp9P3WDST5EkZ+YxrENNRncLVMvRb5OVV8CqwzEs3BNFUmY+bWpW5M2QQFrX9DT84hULFp0WzaAtg7C1smVlr5X4OvsaZ2ApYe9M7QGnbzN4fjW4+hhnbOWB6KWN8B934DOA5NseYnpKKd+71/sNHuBbP4BD82HUMb1uGRuflsv4H06x81wCDau68fFTDWlQVT24u5uc/EJWH7nKgt2XSczIo2WAJ291C6Rd7UqmLs1snU85z5CtQ/Cw92BFrxVUqmDE/1ZnN2uthg6u2iHK6rATs3O3AC/NUsBPgO5CiItA9+LvTScrGcKXQsNn9BreN9JyeO6rgxyKSmZCn3r88HpbFd73UcHOmtD2Ndj7XhcmPV6fqynZDFx8mGErjhKTfP/dFcujup51+TLkSxJzEhm+fThpeWnGGzyoD4RuAytb7bxNdUCExSg7C3l2TIa9n8Lrh8H7Eb1c8mZ6Lv2/OkhyZj7fDGtF4+qqy+Jh5BUUsmx/NF/suIiuUPJKxxq80aW2mn66g4PXD/LGjjcI8gxiYY+FONkacVvYrCRY+xJcPQAd3oYu/wErizgyoMwzxB24+chN07bSDHpcb+GdkJHLCwsPkZSZz4rQliq8S8HexpoRnWqx853OPNbIh3m7LhMyczc/RVxXh1H8QxvfNszoNIMzyWd4a+db5BXmGW9wp0rw8iZo9rI2N772RcjLMN74ygMrGwF+ZBHkpenttJ2kzDwGLDpMfHouy4e0UG1xelLZ1YHP+jdh/Yg2eDjaMerb33lh0SHOxaebujSzEuIXwuR2kzkcf5h3dr+DrkhnvMFt7KDv59BzGlz4VTs4OTXaeOMrD8TyAzw/Cw7Og8Ae4POQ52XeJiUrn4GLDhOXmsOywS0IVq1wehcc4MlPo9oz5YkGnIvPoPecvUzcdJq0bCMGlZnrW6sv41uNJyw2jAn7J1Aki+7/Jn0RAlqP0PYTT4+DRV0her/xxldKzPID/NhyyEmBDqW/+07Nymfg4sNEJ2exZFCwRW5CZSmsrQQvtvZn19udGdDKj68PxdBlZhjfHrlKYZGaVgF44ZEXeLPpm/wc9TNTD001/nRTra4wbCdU8ISVj8PBL9XByWbGsgNclwsHvoCADuDXqlSXSsvW8eKSw1xOzGTxoGDaqpY3o/BwsmPKEw35aVR7ank58f6GUzwxbz/HYu65L1q5MazhMIY0GMK6C+uYdHAShUWFxi2gUm0Y9pv2L9yt78O6l7RnTopZsOwAP7EKMm6Ueu47LUfHS0sPc/FmJgtfak6HQMPtn6LcWX1fN9a92oY5zzchISOXp+cfYOy6EyRk5Jq6NJMSQjCm2RheafgK31/8nnF7x6ErNPJUUwV3bZFPjynakW1fddSW4SsmZ7kBXqiD/bOhajDU6PTQl8nI1TFo6RHO3khn/ovN6FzOl8WbkhCCfk2qsuPtzozoVIufIq7T9dPdLNoTRX6BEeeAzYwQgjebvcmY5mPYEr2Ft3a9RU5BjrGLgLajYMgvUJCvPdw8ukRNqZiY5Qb4qe/g1lXo+O5Db8STlVfAkGVHOR2XxtwBzQgJqqznIpWH4Wxvw7hej7B1dEeCAzyY+stZes3Zw6GoZFOXZlJDGwzlv23+y764fbz222tk5mcavwi/1jBiLwS0h5/HwvfDVKuhCVlmgBcVwt5ZULkh1Hn0oS6RnV/AkOVH+T32Fl+80JRH61fRc5FKadX0cmbZ4BYsfjmYvIIinl94iPfWR3ArO9/UpZnMs3WeZVrHaUQkRBC6LZTUXBM8K3CqBAPXQ9f/wJkNsLCLdnC4YnSWGeBnf4Tki9Dx7Ye6+87JLyR0eTjh0SnM7t+EXg3VBj7mSghBt3qV2TamI692qsn3x+MImbmbH36/Vm4XAfWq0Ys5Xedw+dZlBm8ZzM2sm8YvwspK+9fvy5u0h5qLusKJ1cavo5yzvACXEvbMhIqB2srLB5SrK+SVleEcupLMrOea0LexkXZ+U0rF0c6G93sF8dPI9lTzdGTM2gheWnKE6KTyubdKx2odmd9tPjezbzJoyyBi02NNU0iNjjBiH1QLho2vaVvU5mebppZyyPIC/MJWuHkKOowFqwc7bzJXV8irXx9j/+UkZjzTmCea6n+/cMWw6vm6suG1tvyvX31OxN7i0dl7mLfrUrl8yNmiSguW9FhCli6Ll7e8zMXUi6YpxKUyvLRRuyP//RtY3E071lAxOMsKcClhzwxw99M2oX8AeQWFvL7qOLsvJPLJUw15pnk1AxWpGJq1leDlNgH8NrYTXR/xZsbW8/T5Yi/h0SmmLs3o6leqz/Key7HCisFbBnMq8ZRpCrG20ebEB36vtfYu7Kx2NTQCywrwK7shLhzajQZr2xK/TVdYxMjVv7PzXAJTn2xA/xZ+BixSMZYqbg7Mf7E5i18OJjO3gGcWHOT9DafK3ZL8Wu61WNFrBa52rgzbNowjN46YrpjAblqXSuX68H2odkKWrnz38huSZQX4nk/BxQeaDHygt33wwym2R97kf/3qM7CVv4GKU0ylW73KbB/bidD2NVh79Cohs8rfTofVXKr9eZrPa7+9xq6ru0xXjFs1GPyz1jcevgSW9oCUK6arpwyznAC/ehii92p/KGxLfnjw+mPXWBd+jZFdavNymwDD1aeYlJO9DRP61OPHke3xcXNg1Le/M3jZUWJTys8DNS9HL5Y9uoy6nnUZEzaGzVGbTVeMta22cvP51dpuhl91goi1auGPnllOgO/9FBwrQvPBJX7LxZsZTNh4mtY1PRnTvY7halPMRoOqbvzwelsm9KnH0egUun+2mwW7L6MrLB8POd0d3FnUYxHNKzdn/N7xrDm3xrQFPfIYvLoHvOrCD8O1vVSykkxbUxliGQF+IwIuboPWr4NdyU4oyc4v4PVVx3Gyt+bz55tibaUO1S0vbKytCG1fg9/GdqJ9bS8++fUcfb/Yx/Gr5WODLCdbJ77s9iWdqnVi6uGpLD612LQFeQTA0C3QbZLWRTavlXYOp1JqlhHgB+eBvRu0fKXEb5m46QyXEjP5rH8TvF1LPuWilB2+7hVY9HJzFrzYjFvZOp6ef4D/bDxFWk7Zf8hpb23PrC6zeKzmY8w5PodZx2aZ9pmAlTW0Hw3Dw7ST79cOhB9GQM4t09VUBljGmZg5t+DmaW3/hRJYf+wa73wXwZtdazO2R90HH08pczLzCpi57TwrDkRT0dmeCX3q0beRD+Ih99GxFEWyiI8Of8Ta82vp7t+dKe2m4GjraNqiCvK1duC9M8GlCvSbq+09rtzV3c7EtIwAfwAXb2bw+Nz9NK7uxqphrdXUifI3p66lMf6HU5yKS6NDYCWmPNEA/4pGPDjYBKSUrIxcyaxjs6jlXos5XeZQ3aW6qcuCuGPaXXjSBQgOhR6TSzxFWt6UiwDPzi+g39z9pGbn88ubHdTUiXJHhUWSrw9G8+m2C+gKixjVtTbDO9bCzsYyZhQf1oG4A7y7512EEMzoOIM2vm1MXRLocmDnFG2a1CMAnlyg7Xio/E3ZPpW+mJr3VkrC2kowuJ32kDMkyJtPt12g9+d7OVzGt6ttW7Utax5bg1cFL0b8NoKVZ1aavlfetgI8OhUGbwZZCEt7wvb/qsU/JVRmAnz9sWt8d+wao7rUVifqKCVSxc2BLwc2Z+ngYHLyC+m/8BDvfhdBalbZ3a62umt1vun9DV2qd2FG+Aw+2PcBuQVmEJYB7eG1A9B8EOyfoy3FV6f+3FeZmEJR895KaWXnFzBnx0WW7L2CawVbxvcO4ulmVcvsQ84iWcTCkwuZd2Ie9SvWZ3aX2VRxMpM98S9uh00jITsJOr6nbVz3AFtnlEVldg5czXsr+nQuPp3xG05x/OotWtf0ZMoTDant7Wzqsgxm19VdvL/vfeyt7fms82c0q9zM1CVpslPg1/+DU+vAtyk8sQC8HzF1VSZTZufA1by3ok+PVHFl/Yi2fPRkQyKvp9Nrzh5mbTtPrs7Ip8EbSRe/LqzuvRoXOxdCt4Wy7vw6U5ekcfSEpxfBsysgNUY7SDnsEzU3/g8WHeBq3lsxBCsrwYBWfux4uzOPNfTh852X6Dl7D7vOJZi6NIOo6V6T1Y+tprVPayYfmsykg5Mw+sn3d1P/CXjjsLYkP+xj+LK1NsWiABYc4Lfvc/JWN7XPiaJ/Xi72zH6+Kd+EtsLKSjBk+VGGLj/KlTJ4CpCrnStzu84ltEEo6y+sJ3RbKEk5ZrJnibM3PLtMO77N2hZWPQNrBmqHmpdzFjkHrua9FWPLLyhixYFo5uy4SF5BIaHtazKya22c7W1MXZrebYnewn/3/xcXOxfmdJlDg0oNTF3SXwry4dA82D1d29mw4zvaDqU29qauzKDK1By4mvdWjM3OxopXOtZk5zud6NekKgt2X6brp2Fl8nDlngE9+brX19gIGwb9OogfL/9o6pL+YmMH7cfAG0cgsDvsnAzz28KlHaauzCQsLsDVvLdiSt4uDnz6bGN+eL0tPm4OjFkbwdPzD3DqWpqpS9Orup51WdNnDU28m/DBvg+YdmSa+cyLA7hXh/5fw4vfa3fi3zwF616GtDhTV2ZUFjWFovq9FXNSVCRZf/wa07ecIzkrn/7B1Xnn0bpUci47/5zXFemYGT6TVWdXEeQZxNT2Uwn0CDR1WX9XkAf7P9fODBDW0Ok9betpGztTV6Y3Ft8Hrua9FXOVnqvj898usvxANBXsrBnTrQ4vtfHH1tri/oF7Vzuv7mTSwUlk5GcwsulIBtUbhLWVtanL+rvUGNgyDs7/ApXqwmOfQo2Opq5KLyw+wN/9LoL1x6+xcmhLNXWimKVLCRlM+imSvReTCPR25sPH69OudiVTl6U3KbkpTDk0he0x22ni1YSp7afi52qGB4Sf3wK/vge3YqDBM9rRbq4+pq6qVCw6wNX+3oqlkFKyPfImk3+OJDYlh571q/DBY0FU9zTxHtx6IqXklyu/MPXwVAqKChjTfAz96/bHSpjZvzZ0ObBvNuz7TGs97Pw+tBxusdMqFh3gi/ZEEXYhgZVDW6l5b8Ui5OoKWbw3inm7LlMkJa92rMnwTrXKTNvhzaybTDw4kf1x+7UFQO0mm89eKrdLvqwtyb+0Hdz9oct4aPisdkKQBbHoAAdtD2cV3oqluX4rh49/PcdPEdep6GTHqK61eaGVH/Y2lhUgdyKlZP3F9cw4OgNrYc24luN4vNbj5rcBmJRw6TfYMQniT4FXEIRMgLq9wdxqvYuHDnAhRHVgJVAFKAIWSinnCCE8gbVAABANPCelvOepscY4kUdRzNHxq6lM33KOQ1EpVPOowNjudejXpGqZuCmJzYjlP/v+w/GE43Su3pmJbSZSqYIZzv0XFUHkRtg1FZIvQdVgCPkv1Oxk6sruqzQB7gP4SCmPCyFcgGPAE8BgIEVK+YkQYhzgIaX8v3tdSwW4Up5JKdlzMYnpW85x5no6dSu78O6jdQkJ8ja/u9YHVCSL+Dryaz4//jmOto5MaD2BHgE9TF3WnRUWwIlVsHsapMdBzc5akFdtburK7kpvUyhCiE3A3OKPzlLKG8UhHyalvOcTRhXgiqL1j/986gYzt50nOjmbYH8P3uv5CC1reJq6tFK7fOsy4/eNJzI5kt41ejO+1Xjc7N1MXdad6XLh6GLtcOWcFAjqC10ngJf5NUroJcCFEAHAHqABcFVK6X7br6VKKT3u8J7hwHAAPz+/5jExMQ9cvKKURbrCItaFxzLnt4skZOTRpa4X7/V8hCAfV1OXViq6Ih2LTy1mYcRCPB08+bDth3So1sHUZd1dbjoc+hIOzAVdFjR+ATr9H3j4m7qyP5U6wIUQzsBuYKqUcoMQ4lZJAvx26g5cUf4tJ7+Q5QeimR92iYy8Avo19mVs97r4VbTs1sPI5Eg+2PcBl25d4unApxnTfIz53o0DZCXDvllwZBHIIggeqm2W5ext6spKF+BCCFtgM7BVSjmr+LXzqCkURdGbtGwdC/ZcZtn+KxQWSV5o6cfIrrXxdrHcVcd5hXnMOzGPFWdW4GLnwmuNX+O5us9ha2XGR6SlXdPmx39fBTYO0Po1bcfDCu73f6+BlOYhpgBWoD2wHH3b6zOA5NseYnpKKd+717VUgCvK/d1Mz2XOjousPRqLnbUVoe1rMLxTTVwdzDj07uN8ynlmhM/g8I3DBLgG8E7wO3Ss1tG8H94mXdI6Vs5sAAd3bSFQy1dMckdemgBvD+wFTqG1EQKMBw4D6wA/4CrwrJQy5V7XUgGuKCV3JSmLmdvOs/nkDdwdbRnUJoDBbQPwcLLM1YRSSnZf283M8JlEp0fT2qc177Z4lzoeZn4gy40ICJum7bFibQeNn9fuyCsZb1Mvi1/Ioyjl1em4NGb/doHfziZQwdaaF1r6MaxDDXzdK5i6tIeiK9Sx9vxa5kfMJ1OXyVOBTzGyyUgqVqho6tLuLekiHJwHJ1ZDYZ62EKjtKPBrY/AFQSrAFcXCXbiZwYLdl9l04joCeKJpVUZ0qkltbxdTl/ZQ0vLSmB8xn7Xn1mJvY88rDV/hxXovYm9t5tvxZibC0UXaw86cFK1/vO2bWhuigZboqwBXlDLiWmo2i/deYc3Rq+TqiuhRrzKvda5FU797NoGZrStpV5gZPpPd13ZT1bkqY5qPoYd/D/OeHwfIz4aI1Vr7YeoV8AiA1m9A04Fg56TXoVSAK0oZk5yZx4qDMaw4EE1ajo7WNT15rXNtOgZWMv/wu4OD1w8yI3wGF1Mv0sy7Ge+1eI/6leqbuqz7KyqEcz/Dgc/h2lGo4AEthmkPPfX0wFMFuKKUUVl5BXx75CqL914hPj2Xej6uvNa5Fr0b+ljcXiuFRYVsuLSBub/PJSU3hb41+/JmszfNc6fDO7l6WAvycz8XP/DsD21Glnp1pwpwRSnj8guK2Ph7HAv2XCYqMQv/io4M71iTp5tVw8HWsnY/zMzPZPGpxXwd+TVWwoohDYbwUr2XcLGzkPn+pEtwqPiBZ0Eu1OkJIROhcr2HupwKcEUpJwqLJNsj45kfdpmIa2lUcrZnaPsAnm/hh6eFtSBey7jG7OOz2Rq9FSdbJ56r8xwv1nsRb0fTr44skawkbb+VIwvhxQ3g2+ShLqMCXFHKGSklBy8nM3/3ZfZeTMLO2oreDaswoJU/LQI8LGqe/GzyWZadXsbWmK1YC2v61urL4PqDqeFWw9SllUxBHtg8fHeNCnBFKcfOx2ew+nAMG36PIyO3gEBvZwa08uOpptVwc7ScFZ6xGbGsOLOCjZc2kl+YT5fqXRjacCiNvRqbujSDUgGuKArZ+QVsjrjBqiNXiYi9hb2NFX0a+TKwtR9Nq7tbzF15Sm4Kq8+u5ttz35Ken07zys0Z2mAoHap2sJjfw4NQAa4oyt+cjktj9ZGrbPo9jqz8Qh6p4sLA1v480cQXFwvZdyVbl833F79nZeRK4rPiqe1em6ENhtKzRk/z3jDrAakAVxTljjLzCvjxxHVWHY7hzPV0HO2sebyxLwNb+dOwmhlv/3obXZGOLVe2sPT0Ui7dukQVpyq8XO9lng58Gkdby96WF1SAK4pyH1JKTl5LY/Xhq/wYcZ0cXSENq7oxoJUfjzf2xcnextQl3peUkr1xe1l6einHbh7D1c6VFx55gQFBA/B0sNwTj1SAK4pSYum5Ojb+HseqQ1c5fzMDJztretSvQp9GPnQI9MLOxsrUJd5XRGIES08tZVfsLuys7QjxC6Ff7X60qtIKawPtWWIoKsAVRXlgUkqOX01l7dFYtpyOJz23ALcKtjxavzJ9GvnStlZFbKzNO8yj0qJYfXY1v1z5hYz8DCo7VubxWo/zeK3HCXALMHV5JaICXFGUUskvKGLfpUQ2R9xgW+RNMvMKqOhkR88GVejTyJeWNTzNeul+XmEeYbFhbLq0if3X91Mki2ji1YR+tfvxaMCjZr3KUwW4oih6k6srJOx8IptPXmfH2QRydIV4u9jTu6EPfRv70LS6B1ZmHOYJ2Qn8HPUzGy9tJCotCntre7OeYlEBriiKQWTnF7DzXAI/RVxn1/lE8guK8HVzoE9jX/o08qFhVTez7c2WUnIm+QwbL2006ykWFeCKohhcRq6O387e5KeIG+y9mIiuUOJf0ZHHGvoQEuRN42ruZjtnbs5TLCrAFUUxqrRsHVvPxPPTyescuJxMYZHErYIt7QMr0amOF53reOHt6mDqMu/oTlMsrXxa0alaJzpV60Rlp8pGrUcFuKIoJpOWrWPfpSR2X0gg7HwiCRl5AAT5uNK5rhed6njR3N8DWzO7O/9jimVz1GbCYsOIy4wDIMgziM7VO9OpeieCPIOwEoatWwW4oihmQUrJufgMws4nsvtCAuHRqRQUSZztbWhXuyKd6njTua6X2R3aLKXk8q3LhF0LY3fsbiISI5BIvCt406FaBzpX70wrn1ZUsNF/3SrAFUUxSxm5Og5cTtYC/XwC19NyAQj0di6+O/emRQ0P7G3MqzMkJTeFfXH7CIsN48D1A2TpsrC3tqe1T2s6Ve9Ex6od9TbVogJcURSzJ6XkUkImuy8kEnY+kSNXUsgvLKKCrTVNqrsTHOBBcIAnTf3ccTWjDbd0hTqO3jzK7tjd7L62+45TLfU86z10N44KcEVRLE52fgEHLyez92IS4TEpRF5Pp0iCEPBIFVeC/T3+DPWqZjLlIqXk0q1L7L62+29TLTM7zaRHQI+HuqYKcEVRLF5mXgEnrt4iPCaFYzGpHI9JJSu/EAAfNwea+3vQIsCT5v4eBPm4msXK0JTcFPZe20sXvy642rk+1DVUgCuKUuYUFBZxLj6DYzGpHI1OITw6lfh0bQ7dyc6aZv4eNPf3INjfkwZVXXF3tKwzQf+gAlxRlDJPSkncrRyOxaQSHq2F+vmbGfwRc75uDtTzdaWerxv1fFyp7+tKNY8KZrtS9A93C3Dz3+BXURSlhIQQVPNwpJqHI/2aVAW0rXEjYm8ReT2dyBvpRF5PZ+e5BIqKQ93F3oYgX1fq+bhq4e7jSmBlZ7PrerkTFeCKopRprg62dAj0okOg15+v5eQXcv5mRnGopxF5PZ21R2PJ0Wnz6TZWgtrezn8Gej0fV2p5O+PtYm9Wd+sqwBVFKXcq2GltiU2qu//5WmGRJCY568+79Mgb6ey7mMSG43F//oyTnTU1vJyoWcmZGpWcqPnH115OOJvgxCIV4IqiKIC1laCmlzM1vZzp08j3z9cTM/I4F5/OlaQsohKziErK4vjVVH46eZ3bHyF6u9j/PdSLv67u6WiwLQJUgCuKotyDl4s9Xi5/n4IBbU/0qynZRCVmElUc7leSsth65iYpWbF//py1lcDP05GPnmxIm1oV9VqbCnBFUZSH4GBrTZ3KLtSp/O9tZm9l5xOVlMWVxCyikjK5kpRFRWf9tzCqAFcURdEzd0c7mvnZ0czPw6DjmNfejYqiKEqJqQBXFEWxUCrAFUVRLFSpAlwI0VMIcV4IcUkIMU5fRSmKoij399ABLoSwBuYBvYB6wAtCiHr6KkxRFEW5t9LcgbcELkkpo6SU+cAaoJ9+ylIURVHupzQBXhWIve37a8Wv/Y0QYrgQIlwIEZ6YmFiK4RRFUZTblSbA77Sjy7/2ppVSLpRSBkspg728vO7wFkVRFOVhlGYhzzWg+m3fVwOu3+sNx44dSxJCxDzkeJWApId8r6mp2k3DUmu31LpB1W4o/nd68aEPdBBC2AAXgBAgDjgKDJBSnnnYCu8zXvidNjS3BKp207DU2i21blC1G9tD34FLKQuEECOBrYA1sNRQ4a0oiqL8W6n2QpFS/gL8oqdaFEVRlAdgSSsxF5q6gFJQtZuGpdZuqXWDqt2ojHqosaIoiqI/lnQHriiKotxGBbiiKIqFMosAF0JUF0LsEkKcFUKcEUK8Vfy6pxBiuxDiYvFnj9ve837xJlrnhRCPmq76P+uxFkL8LoTYXPy9RdQuhHAXQqwXQpwr/u/fxoJqH1P85+W0EOJbIYSDudYuhFgqhEgQQpy+7bUHrlUI0VwIcar41z4XRjgi/S61zyj+M3NSCPGDEML9tl8zi9rvVPdtv/aOEEIKISqZW90PREpp8g/AB2hW/LULb1NTqQAAA7dJREFUWn95PWA6MK749XHAtOKv6wERgD1QA7gMWJv49zAWWA1sLv7eImoHVgDDir+2A9wtoXa0bRuuABWKv18HDDbX2oGOQDPg9G2vPXCtwBGgDdpK6F+BXiaqvQdgU/z1NHOs/U51F79eHa39OQaoZG51P8iHWdyBSylvSCmPF3+dAZxF+x+0H1rAUPz5ieKv+wFrpJR5UsorwCW0zbVMQghRDXgMWHzby2ZfuxDCFe0P+RIAKWW+lPIWFlB7MRugQvGiMke0lcBmWbuUcg+Q8o+XH6hWIYQP4CqlPCi1ZFl523uMWruUcpuUsqD420NoK7HNqva7/DcH+Ax4j79v/WE2dT8Iswjw2wkhAoCmwGGgspTyBmghD3gX/1iJNtIyotlofyCKbnvNEmqvCSQCy4qnfxYLIZywgNqllHHAp8BV4AaQJqXchgXUfpsHrbVq8df/fN3UhqLdmYKZ1y6EeByIk1JG/OOXzLruuzGrABdCOAPfA6OllOn3+tE7vGaSfkghRB8gQUp5rKRvucNrpurltEH7J+Z8KWVTIAvtn/J3Yza1F88X90P75+7/t3P2rFFEURh+3kKEWKkgIikiIrYWFkEtxFiISOqAwS38ESILgn/AysLGSoONBLXXXhDRKH6gYsAtopaCTYRjcc+SQTYLK5i9F94Hhtk9dy48d9k5w5x7Zw4BeyQtj+syIlbrGtrtXKsbg6Q+8BtYGYZGHFaFu6QZoA9cH9U8IlaF9ziqSeCSdlGS90pErGb4W97CkPvvGZ/4RVr/kVPAoqR1yjvRz0q6RxvuA2AQEc/y+wNKQm/B/RzwJSJ+RMQmsAqcpA33IZO6DtgqVXTjU0FSD7gIXMryAtTtfoRywX+V5+ss8ELSQer23pYqEnjO6t4B3kXEzU7TY6CXn3vAo058SdJuSYeBo5SJhh0nIq5FxGxEzAFLwNOIWKYN9w3gq6RjGVoA3tKAO6V0Mi9pJv8/C5S5kxbch0zkmmWWn5Lmc8yXO312FEnngavAYkT86jRV6x4RryPiQETM5fk6oCye2KjZeyzTnkXNC/dpym3JGvAytwvAfuAJ8DH3+zp9+pSZ4g9UMisMnGFrFUoT7sBx4Hn+9g+BvQ253wDeA2+Au5QVBFW6A/cptfpNSuK48i+uwIkc72fgFvk09RTcP1FqxsPz9XZt7qO8/2pfJ1eh1OQ9yeZH6Y0xplGqKKEYY4yZHCdwY4xpFCdwY4xpFCdwY4xpFCdwY4xpFCdwY4xpFCdwY4xplD96cklJaEu9GAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The reference state\n",
    "X = np.array([250, .5])\n",
    "plt.plot(nu,greybody(nu,X))\n",
    "\n",
    "# A different temperature\n",
    "X = np.array([300, .5])\n",
    "plt.plot(nu,greybody(nu,X))\n",
    "\n",
    "# A different emissivity\n",
    "X = np.array([250, 1])\n",
    "plt.plot(nu,greybody(nu,X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Simulate an \"observed spectrum\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x123707c88>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZd7/8fd3Jo00SCdAINI7IYSOWBALUlQWaSIoLo9lV1x3V3Fdn9W1/NRnlRUVlUUlIooIKsWKqFhAIQFCJ5GaRgg9hfT790dGl3VBEsjkzJl8X9eV68ycyWQ+A8knJ/cptxhjUEopZT8OqwMopZQ6P1rgSillU1rgSillU1rgSillU1rgSillUz71+WKRkZEmPj6+Pl9SKaVsLzU19bAxJuqX6+u1wOPj40lJSanPl1RKKdsTkf1nWq9DKEopZVNa4EopZVNa4EopZVNa4EopZVNa4EopZVNa4EopZVNa4EopZVP1ehy48nyVVYajRWUcLizlSGH18nBhKSLCtd1iado4wOqISikXLfAGwBhD7okS8gtKfy7m/MLS/yjpn5ZHi8s42yXiH/9wO4PbRzE2KY4hnWLw89E/4JSykha4l8s9cYr7Fm/mm4zD//VYkJ+TyBB/IoP9aRURSK/4MCKD/IgM8SciyJ/IYD8igv2JCvbnWHEZi1OzWJyaxR0LNhAR5Mf1PZsztncc7WJCLHhnSimpzxl5kpKSjJ5KXz+MMSzdlMNDS7dSUWn43eVt6RAT4ipnPyKD/Wnk56z1162sMnydns876zP5fEceFVWGhLgmjO0dx/DusYQE+Lrh3SjVsIlIqjEm6b/Wa4F7n6NFZfz1gy18tOUgvVqF8cyYHsRHBtX56xwpLOX9jdm8sz6TjEOFNPJ1cm33WG5MiqN3fBgiUuevqVRDpAXeQHyxM4/7l2zheHEZ9w7twLTBrXE63Fukxhg2ZR5nUUomy9NyKSytoHVkEGOS4hjdqznRIbrjU6kLoQXu5QpLK3hsxXYWrs+kY9MQnr0xgc7NQus9R3FZBR9tOcii9Zms23cUp0O4rEM0E/u15LIO0fWeRylvoAXuxdbtPcof391E9rFTTBvchj8MbYe/T+3Ht+vanvxC3k3NYklqFocKSvlNrxb8fVQXAv1037lStaEF7oVKyiuZuTKdOd/sIS4skGdv7EFSfLjVsf5LRWUVz3/xI7O+yKBNVDAvTkikQ1M9ckWpmjpbgeuBvDa1LecEo174jle+3sP4Pi35ePrFHlneAD5OB38Y2p4FU/tyvLicUS9+yzvrD1CfGw9KeSMtcJupqKzixS9/5LoXv+NYcRmv39KbJ67vRpC/5w9LDGgbyUfTB9GrVRj3L9nCH97ZRFFphdWxlLItz/+pVz/be7iIexdtYuOB4wzvHsujo7oSFuRndaxaiQ4J4I1b+/Lilz/yz8/T2Zx1ghcmJFqyw1Upu9MtcBswxjB/7T6GPfcNe/KLmDW+Jy9MSLRdef/E6RDuHtKOt37bj6KyCq6b/R0LftivQypK1ZIWuA08umIHDy3dRu+Lwvn0nsGM7NHM6kh1ol/rCD66+2L6tY7gwfe38vu3N1JQUm51LKVsQwvcwy1Py+G17/YyZUA8ybf09rqrAUYE+zNvSm/uu7oDH289yIjnv2Vr9gmrYyllC1rgHuzHQwXcv2QzSa3CePDaTl57arrDIdx5aVsWTutHaUUVN8xeQ/KafTqkotQ5aIF7qKLSCm5/cwOBfk5emJCIr9P7/6t6x4fz4d0XM6hdJH9bto07F2zgxCkdUlHqbLy/FWzIGMMD721hT34hs8b19Lphk18THuTH3JuT+Muwjqzcnsfw578hLfO41bGU8kha4B5o/vf7WZaWwx+v7MCAtpFWx6l3DocwbXAbFt3en6oq+M3La3jt2706pKLUL5yzwEWkg4hsOu3jpIjcIyLhIrJSRDJcy7D6COztNh44xqMrtjOkYzR3XNLG6jiWSmwZxod3D+LSDtH8fcV2Zn6eYXUkpTzKOQvcGLPLGJNgjEkAegHFwPvADGCVMaYdsMp1X12Ao0Vl3LVgAzGhATx7YwION18G1g6aBPrxyk29uDGpBbNWZTBrlZa4Uj+p7ZmYQ4Ddxpj9IjIKuNS1Phn4Cri/7qI1LJVVhnve2cThwjKW3DGAxoE6s81PHA7hyRu6U1FleHZlOk6HcNdlba2OpZTlalvg44C3XbdjjDG5AMaYXBHRiz1fgOe/yODr9HyeuL4b3Vo0tjqOx3E4hP/7TQ8qqwz/9+kufJ3V4+RKNWQ1LnAR8QNGAg/U5gVEZBowDaBly5a1CtdQrE7P57lVGdyQ2JzxfeKsjuOxnA7hmTHVJf7ERztxOhxMHXSR1bGUskxttsCvATYYY/Jc9/NEJNa19R0LHDrTk4wxc4A5UH098AtK64Wyj5/inoUb6RATwuPXdfPak3Xqio/TwcyxCVRWGR5dsR0fhzB5QLzVsZSyRG0OIxzPv4dPAJYBk123JwNL6ypUQ1FWUcVdCzZQXmmYPTHxvGaJb4h8nQ5mje/J0M4x/G3ZNt78fr/VkZSyRI0KXEQCgaHAe6etfhIYKiIZrseerPt43u3xD7ezKfM4/xjTndZRwVbHsRVfp4MXJyQypGM0f/1gKwvXHbA6klL1rkZDKMaYYiDiF+uOUH1UijoPSzdlk7x2P7cNuoiru8ZaHceW/HwczL4pkWlvpPLA+1twOoQxSboPQTUceiamBTLyCpixZAu948O4/5qOVsexNX8fJ69M6sWgtpHct2Qz72/MsjqSUvVGC7yeFZZWcPubqQT5+zSYi1S5W4CvkzmTkuh3UQR/XJTG0k3ZVkdSql5oe9QjYwwzlmxm7+Einh/fk5jQhnORKndr5Ofk1SlJJMWHc++iND7cnGt1JKXcTgu8HiWv2ceKzbn86aoO9G8Tce4nqFoJ9PPh9Sm96RnXhLsXbuSTrQetjqSUW2mB15MNB47x+Ec7uKJTNLfrGYRuE+Tvw+u39KZ7i8b87q0NrNyed+4nKWVTWuD14EhhKXct2EDTxgE8M0YvUuVuIQG+JN/ahy7NQrlzQSpf7NQSV95JC7wePPHRTo4UlvHSxF56kap6Ehrgyxu39qVD0xBun7+B1en5VkdSqs5pgbvZj4cKeX9jFjf3b0XX5nqRqvrUONCXN6f2pU10MNPeSGHDgWNWR1KqTmmBu9lzqzII8HVy+6U67m2FJoF+vDm1DzGhAUx7I4XMo8VWR1KqzmiBu9GugwWs2JzD5AHxRAb7Wx2nwYoI9ue1Kb0prahiavJ6TpboRMnKO2iBu9HMlekE+fkw7eLWVkdp8NpGB/PyTb3Yk1/EXQs2UFFZZXUkpS6YFribbM0+wSfbDnLroIsIC/KzOo4CBraN5LHruvJNxmH+tmybTpKsbK+2M/KoGpq5Mp3QAB+dcMDDjOvTkr1Hinhl9R5aRwXr/4+yNd0Cd4ONB46xauchpg1uTeNGetigp7n/qo5c3aUpj324nc/1RB9lY1rgbvDsynTCAn2ZMlC37jyRwyHMHJtAt+aNuXvhRrZmn7A6klLnRQu8jq3fd5RvMg5z+yVtCPbXESpP1cjPydybk2jcyJfbklM4eKLE6khK1ZoWeB179rN0IoP9ubl/vNVR1DlEhwbw6uTeFJSUMzV5PUWlFVZHUqpWtMDr0Jrdh1m75wh3XtpG57e0ic7NQnl+Qk925J7knnc2UVmlR6Yo+9ACryPGGJ79LJ2moQFM6NvS6jiqFi7vGMP/Du/Myu15PPXJTqvjKFVjOkhbR77OOEzK/mM8el1XAnx169tupgy8iL2Hi5jz9R7iI4L0l7CyBS3wOlC99b2L5k0aMVYn1bWth4Z3Zv/RYh5aupWW4YEMahdpdSSlfpUOodSBVTsOkZZ1gruHtMXPR/9J7crH6eD58T1pFx3MHQtSycgrsDqSUr9K2+YCVVUZnl2ZTquIQG5IbGF1HHWBQgJ8eXVKb/x9nNyavJ7DhaVWR1LqrLTAL9Cn2w6yPfck04e00xnmvUTzJo2YOzmJQydLmfZGCiXllVZHUuqMatQ4ItJERBaLyE4R2SEi/UUkXERWikiGaxnm7rCeprLKMPPzdNpEBTEqobnVcVQdSohrwsyxCWw4cJw/L96sF75SHqmmm4zPAZ8YYzoCPYAdwAxglTGmHbDKdb9BWbE5h/S8Qu65oj1OnefS6wzrFst9V3dgeVoOMz/PsDqOUv/lnAUuIqHAYOBVAGNMmTHmODAKSHZ9WjJwnbtCeqKKyiqe+zyDjk1DuLZbrNVxlJvccUkbbkxqwaxVGSzdlG11HKX+Q022wFsD+cDrIrJRROaKSBAQY4zJBXAto8/0ZBGZJiIpIpKSn+89E8t+sCmHPYeLuOeK9jrLvBcTER67rht94sO5b/Fm0jKPWx1JqZ/VpMB9gETgJWNMT6CIWgyXGGPmGGOSjDFJUVFR5xnTs5RXVjFrVQZdm4dyVZcYq+MoN/PzcfDSTYlEBvszbX4KeSf1wlfKM9SkwLOALGPMD677i6ku9DwRiQVwLQ+5J6LnWZyaxYGjxdw7tD0iuvXdEEQE+zN3chIFJRVMm5+qR6Yoj3DOAjfGHAQyRaSDa9UQYDuwDJjsWjcZWOqWhB6mtKKS51dlkBDXhMs6nHHUSHmpTrGhPHtjAmmZx5mxRI9MUdar6an0vwcWiIgfsAe4heryXyQiU4EDwBj3RPQsi9ZnknOihCdHd9et7wbo6q5N+dOV7fnHZ+l0aBrKHZe2sTqSasBqVODGmE1A0hkeGlK3cTxbSXklL3z5I73jw7hYr5PRYN11WVt25RXy9Kc7aRcdzBWddT+IsoaeOlgLC344QN7JUu4d2kG3vhswEeHp0d3p2qwx0xduJF2vmaIsogVeQ8VlFbz01Y8MaBNB/zYRVsdRFmvk5+RfNycR6O/DbckpHC0qszqSaoC0wGvojbX7OVxYxh+vbG91FOUhmjYOYM6kXhw8WcKdC1Ipr6yyOpJqYLTAa6CwtIJXVu/mkvZR9GoVbnUc5UF6tgzjqdHd+H7PUR5ets3qOKqB0QkdauCNtfs4VlzOvUN161v9t+t7tmDnwQJeWb2Hjk1DmKQTWqt6olvg51BVZVjw/QEGtImgR1wTq+MoD3XfVR25vGM0Dy/fzpofD1sdRzUQWuDn8M2Ph8k+forxfXSORHV2Tofw3LgEWkcGcedbG9h/pMjqSKoB0AI/h4XrDhAW6MuVes0TdQ4hAb7MnVx9usTU5BQKSsotTqS8nRb4r8gvKGXl9jxGJ7bA30dnmlfn1ioiiNkTEtl7uIjpCzdRWaWn2yv30QL/FUs2ZFFRZRjXR2eaVzU3oG0kD4/ozBc7D/H0pzutjqO8mB6FchbGGBauO0Cf+HDaRodYHUfZzKT+8ezKqz4ypUNMiE54rdxCt8DPYu2eI+w7Uqxb3+q8/W1EF/q3jmDGe1vYeOCY1XGUF9ICP4uF6zIJDfBhmE6Xps6Tr9PB7ImJNA0NYNr8VJ0IQtU5LfAzOFZUxidbD3J9z+YE+OrOS3X+woL8mDs5iaLSCu5csIGyCj3dXtUdLfAzeG9jNmWVVYzvq8d+qwvXPiaEp0Z3J3X/MZ74aIfVcZQX0QL/hZ92XibENaFj01Cr4ygvMaJHM6YOuoh5a/bp7PaqzmiB/0Lq/mNkHCpkvO68VHVsxjUd6RMfzowlW9h58KTVcZQX0AL/hbfXZRLk52R492ZWR1Fextfp4IUJPQkO8OH2+amc1DM11QXSAj/NiVPlfLglh5EJzQny10PkVd2LDg1g9sREso6d4o+L0qjSMzXVBdACP82yTdmUlFcxQS9cpdyod3w4fxnWiZXb83hp9W6r4ygb0wJ3Mcbw1rpMujQLpVuLxlbHUV7uloHxjOjRjGc+28W3GXr5WXV+tMBdNmedYEfuScbp1reqByLCkzd0o210MHcv3Ej28VNWR1I2pAXusnD9ARr5OhmVoDsvVf0I8vfh5Zt6UVZRxZ1vplJaUWl1JGUzNSpwEdknIltEZJOIpLjWhYvIShHJcC3D3BvVfYpKK1i2KYdru8cSGuBrdRzVgLSOCuYfY3qQlnWCR5ZvtzqOspnabIFfZoxJMMYkue7PAFYZY9oBq1z3bWl5Wg5FZZU6646yxNVdm3LHpW1464cDLErJtDqOspELGUIZBSS7bicD1114HGu8vT6T9jHBJLbUOS+VNf44tD0D20bw1w+2sjX7hNVxlE3UtMAN8JmIpIrINNe6GGNMLoBrGX2mJ4rINBFJEZGU/Pz8C09cx7bnnCQt8zjjerdERKyOoxooH6eDWeN6EhHkx+1vpnK8uMzqSMoGalrgA40xicA1wF0iMrimL2CMmWOMSTLGJEVFRZ1XSHdauP4Afj4ObkhsbnUU1cBFBPsze2IieSdLmL5wk57ko86pRgVujMlxLQ8B7wN9gDwRiQVwLQ+5K6S7nCqr5P2N2Qzr2pQmgX5Wx1GKni3D+NuILqxOz+e5VRlWx1Ee7pwFLiJBIhLy023gSmArsAyY7Pq0ycBSd4V0l4+25FJQUqHHfiuPMrFvS0YntuC5VRl8udN220WqHtVkCzwG+FZE0oB1wIfGmE+AJ4GhIpIBDHXdt5WF6w/QOjKIvheFWx1FqZ+JCI9f35XOsaFMX7iRA0eKrY6kPNQ5C9wYs8cY08P10cUY87hr/RFjzBBjTDvX8qj749adjLwC1u87xtjecbrzUnmcAF8nL9/UC4Db30ylpFxP8lH/rcGeiblwfSa+TmF0L50tXHmmlhGBPDeuJ9tzT/Lg+1sxRndqqv/UIAu8tKKS9zZkcWXnpkQG+1sdR6mzuqxjNNOHtGPJhizmrdlndRzlYRpkgX+6LY9jxeWM01l3lA1MH9KOoZ1jeHTFdr5O97xzKZR1GmSBv/3DAeLCGzGwTaTVUZQ6J4dD+OfYBNrHhHDXWxvYnV9odSTlIRpcge87XMTaPUcYmxSHw6E7L5U9BPn7MHdyEn5OB7clp3CiWKdjUw2wwBeuz8TpEMYk6fCJspcWYYG8MqkXWceKueutDZRXVlkdSVmsQRV4eWUVi1OzuLxjNDGhAVbHUarWkuLDeeL6bnz742EeW6GXn23oGtTMvat25HG4sJTxuvNS2diYpDjS8wr41zd7aRcTwk39WlkdSVmkQW2Bv7Uuk9jGAVzS/owXTlTKNmZc04nLOkTx8LJtrNmtc2o2VA2mwDOPFvNNRj5jkuJw6s5LZXNOhzBrfE8uigzizgUb2H+kyOpIygINpsDfdc10Mra3Dp8o7xAS4MvcydUTZE1NTuFkiR6Z0tA0iAKvqKxiUUoWl7SPonmTRlbHUarOtIoIYvbERPYdLuLutzdSqdcQb1AaRIGvTs/n4MkSxvXWy8Yq7zOgTSSPjOrCV7vyefLjHVbHUfWoQRyFsjg1i8hgP4Z00p2XyjtN7NuK9IP/PjLlRj3PoUHw+i3wgpJyvth5iOHdm+Hr9Pq3qxqwh4Z3ZlDbSB58fwsp+2x1dWd1nry+0VZuz6O0oooRPWKtjqKUW/k4Hbw4IZEWYYH8z/xUMo/qRBDezusLfHlaDs2bNCKxZZjVUZRyu8aB1UemlFVW8ds3UigqrbA6knIjry7wY0VlfJNxmOE9YnXWHdVgtIkK5sUJiaTnFXDPOzq7vTfz6gL/eOtBKqoMI3s0szqKUvVqcPsoHhremZXb83hm5S6r4yg38eqjUJalZdM6KojOsaFWR1Gq3k0ZEE96XiEvfrmb9jEhjEpobnUkVce8dgs872QJP+w9ysgezXT4RDVIIsIjI7vQ96Jw/rx4M6n79cgUb+O1Bb5icy7GwPDuOnyiGi4/Hwcv3dSL5k0acVtyCnt0Nh+v4rUFvjwth86xobSNDrY6ilKWCg/yY94tvXGIMOX19RwuLLU6kqojNS5wEXGKyEYRWeG6Hy4iK0Ukw7X0mOP0Mo8WsynzOCMTdOtbKai+ZsqrU3pzqKCEqfPWU1ymhxd6g9psgU8HTr/QwgxglTGmHbDKdd8jLEvLAWB4dz15R6mfJMQ14fnxiWzJPqEXvvISNSpwEWkBXAvMPW31KCDZdTsZuK5uo52/5Wk59GoVRouwQKujKOVRhnaO4ZGRXfh8xyEeXrYNY7TE7aymW+D/BO4DTp9FNcYYkwvgWp7xSlEiMk1EUkQkJT8//4LC1kRGXgE7DxYwQre+lTqjSf3j+Z9LWjP/+/288vUeq+OoC3DOAheR4cAhY0zq+byAMWaOMSbJGJMUFRV1Pl+iVpan5eAQGKYFrtRZ3X9VR0b0aMaTH+9k6aZsq+Oo81STE3kGAiNFZBgQAISKyJtAnojEGmNyRSQWOOTOoDVhjGFZWg7920QQHaKzzit1Ng6H8I8x3Tl0soQ/v7uZmNAA+rWOsDqWqqVzboEbYx4wxrQwxsQD44AvjDE3AcuAya5PmwwsdVvKGtqafZJ9R4r11HmlasDfx8mcSUm0jAhk2hspZOQVWB1J1dKFHAf+JDBURDKAoa77llqWlo2vU7i6iw6fKFUTjQN9mXdLb/x9nUx5fT15J0usjqRqoVYFboz5yhgz3HX7iDFmiDGmnWtp6Xm6VVWGFZtzuaR9FI0Dfa2MopSttAgL5PUpvTlWXMYtr6+nUC9BaxtecyZmyv5j5J4oYYQOnyhVa12bN2b2xER25RVw54INlFdWnftJynJeU+DL03II8HVwRacYq6MoZUuXdojmieu78nV6Pg++v0WPEbcBr7icbEVlFR9tyWVIpxiC/L3iLSllibG9W5J97BSzvviR5k0CmX5FO6sjqV/hFW23ZvcRjhSVMUKvPKjUBfvD0PZkHy9h5ufpNGsSwBid4d5jeUWBL0vLIcTfh0s7uP9EIaW8nYjw/27oRt7JEh54bwsxoQEMbq8/W57I9mPgpRWVfLr1IFd2aUqAr9PqOEp5herriCfSNjqYOxdsYHvOSasjqTOwfYF/tSufgtIKvXSsUnUsJMCXebf0ISTAh1vmrSPzaLHVkdQv2L7Al6flEB7kx4A2ehqwUnWtaeMA5t3Sh5LyKibO/YGDJ/REH09i6wIvLqtg1Y5DDOvWFF+nrd+KUh6rQ9MQkm/tw9GiMibO/V5n9PEgtm69ldvzOFVeqUefKOVmCXFNeG1Kb7KPn2LSq+s4XlxmdSSFzQt8eVouTUMD6B0fbnUUpbxen4vC+dfNSew+VMjk19dTUFJudaQGz7YFfqK4nNXphxjePRaHQ6yOo1SDcHG7KF6cmMi27BNMnZfCqbJKqyM1aLYt8E+3HaS80ujRJ0rVs6GdY5g5NoGU/UeZNj+F0gotcavYtsCXpeXQKiKQbs0bWx1FqQZnRI9mPDW6O99kHOauBRv14lcWsWWB5xeUsmb3YUZ0b4aIDp8oZYUxSXE8OqoLn+/I4w/vbNJZ7i1gy1PpP9qSS5VBh0+Ustik/vGcKq/kiY92EuDr5OnR3XWfVD2yZYEvT8uhQ0wI7WNCrI6iVIM3bXAbissq+efnGTTydfL3UV30L+N6YrsCzz5+ipT9x/jzVR2sjqKUcpk+pB2nyip55es9BPo5mXFNRy3xemC7Al+RlgPA8O4676VSnkJEmHFNR06VV5d4Iz8n91zR3upYXs92Bb4sLYcecU1oFRFkdRSl1GlEhIdHdPl5OCXQz8m0wW2sjuXVbFXgu/ML2ZZzkr9e28nqKEqpM3A4hKdGd6fEtWOzka+TSf3jrY7ltWxV4MvTchCB4XrtE6U8ltMhzBybQEl5FQ8t3UaAr1Nn9XET2xwHboxheVoOfeLDado4wOo4Sqlf4et08MKEnlzcLpL7l2xmuWvflapb5yxwEQkQkXUikiYi20TkEdf6cBFZKSIZrmWYO4Nuzz3J7vwiPfZbKZsI8HUyZ1ISSa3C+cM7m1ixWUu8rtVkC7wUuNwY0wNIAK4WkX7ADGCVMaYdsMp1322Wp+XidAjXdNWjT5Syi0Z+Tl6dkkTPlk24++2NLE7NsjqSVzlngZtqha67vq4PA4wCkl3rk4Hr3JKQfw+fDGobSXiQn7teRinlBiEBviTf2ocBbSL507tpzF+7z+pIXqNGY+Ai4hSRTcAhYKUx5gcgxhiTC+BaRp/ludNEJEVEUvLz888r5IYDx8k+foqRPXT4RCk7CvTzYe7kJK7oFM1DS7cx5+vdVkfyCjUqcGNMpTEmAWgB9BGRrjV9AWPMHGNMkjEmKSoq6rxCLk/Lwc/HwZVdYs7r+Uop6wX4Onnppl4M7x7LEx/tZObKdIzRC2BdiFodRmiMOS4iXwFXA3kiEmuMyRWRWKq3zt3iup7NaR8TQkiAr7teQilVD3ydDp4b15MAXyfPrcrgVHklD+hp9+ftnAUuIlFAuau8GwFXAE8By4DJwJOu5VJ3hUyIa0JCXBN3fXmlVD1yOoSnR3cn0M/JnK/3UFxWwd9HdtWrGJ6HmmyBxwLJIuKkeshlkTFmhYisBRaJyFTgADDGjTmVUl7E4RAeGdmFRn5OXlm9h+KySp4e3R0fp21OTfEI5yxwY8xmoOcZ1h8BhrgjlFLK+4kIM67uSJCfD8+uTKe0vIqZYxPw89ESrylbnUqvlPIuIsLdQ9oR6OfksQ93cKq8ktkTEwnwdVodzRb0V51SynK3Xdyax6/vype7DnHrvPUUlVZYHckWtMCVUh5hYt9WPDOmB9/vOcLNr63jZEm51ZE8nha4Uspj3JDYghcnJLI56zgT/vU9R4vKrI7k0bTAlVIe5ZpuscyZlER6XiHj5qzlUEGJ1ZE8lha4UsrjXNYxmnlTepN17BQ3vryW7OOnrI7kkbTAlVIeaUDbSOZP7cORojJufHktGXkFVkfyOFrgSimP1atVOG//th+lFVXcMHsNX+1y2xU7bEkLXCnl0bo2b8zS3w2kRXggt85bz7zv9upFsFy0wJVSHq95k0Ysvr0/l3eM4eHl2/nrB1spr6yyOpbltMCVUrYQ5O/DnEm9uP2SNiz44QCTX1vH8eKGfZihFrhSyjYcDmHGNR35x5gepOw7xvWz17A7v/DcT/RSWuBKKdv5Ta8WvPXbvpw8Vc71L37HtxmHrY5kCS1wpZQtJcWH88FdA4lt3HdjUBAAAAk8SURBVIjJr69j/vf7rY5U77TAlVK2FRceyOI7+nNJ+yge+mArf1u6lYoGtHNTC1wpZWshAb786+Ykpg1uTfLa/dwybz0nTjWMC2FpgSulbM/pEP4yrBNPj+7O93uOcP3s79h3uMjqWG6nBa6U8ho39o5j/tS+HCsqY9SL37Fmt3fv3NQCV0p5lX6tI1h61yCiQ/y5+dV1vL3ugNWR3EYLXCnldVpGBLLkzgEMbBvJA+9t4ZHl27xy56YWuFLKK4UG+PLq5CRuHXgRr3+3j7FzvifzaLHVseqUFrhSymv5OB3874jOzBrfk/SDBQyb9Q3L03KsjlVntMCVUl5vZI9mfDT9YtpGB/P7tzfy53fTvGLi5HMWuIjEiciXIrJDRLaJyHTX+nARWSkiGa5lmPvjKqXU+YkLD2TR//Tnd5e1ZfGGLEY8/y1bs09YHeuC1GQLvAL4ozGmE9APuEtEOgMzgFXGmHbAKtd9pZTyWL5OB3+6qgNv3daP4rJKrp/9HXO/2UNVlT2vL37OAjfG5BpjNrhuFwA7gObAKCDZ9WnJwHXuCqmUUnWpf5sIPp5+MZd2iOaxD3dwy7z15BeUWh2r1mo1Bi4i8UBP4AcgxhiTC9UlD0Sf5TnTRCRFRFLy8/MvLK1SStWRsCA/5kzqxaPXdeX7PUe45rmvWZ1ur46qcYGLSDCwBLjHGHOyps8zxswxxiQZY5KioqLOJ6NSSrmFiDCpXyuW/W4QEUH+TH5tHY+t2E5pRaXV0WqkRgUuIr5Ul/cCY8x7rtV5IhLrejwW0NlGlVK21KFpCEt/N5Cb+7di7rd7ucEmE0XU5CgUAV4Fdhhjnj3toWXAZNftycDSuo+nlFL1I8DXyd9HdeVfNyeRc/wUw2d9y6L1mR49gXJNtsAHApOAy0Vkk+tjGPAkMFREMoChrvtKKWVrQzvH8PH0wSTENeG+JZv5/dsbPfbytFKfv12SkpJMSkpKvb2eUkqdr8oqw8urd/PsynSahgbw5OhuXNzOmv14IpJqjEn65Xo9E1Mppc7A6RDuuqwti2/vj5+Pg0mvruPutzdyqKDE6mg/0wJXSqlf0bNlGB9Pv5h7rmjHJ1sPMuSZ1cxfu49KDzj5RwtcKaXOIcDXyT1XtOeTey6me4vGPLR0Gze8tMbyU/G1wJVSqoZaRwXz5tS+PDcugexjxYx84VseXbGdQosujKUFrpRStSAijEpozqp7L2V8n5a89t1ernhmNZ9sza33Qw61wJVS6jw0DvTl8eu7seSOAYQF+XH7mxu4LTmlXieN0AJXSqkLkNgyjOW/G8iDwzqxds8Rhs5czUtf7aa8HqZw0wJXSqkL5ON08NvBrVl57yUMbhfFU5/sZPisb0nZd9Str6sFrpRSdaR5k0bMuTmJf92cREFJOb95eS0zlmzmWFGZW15PC1wpperY0M4xrLz3EqYNbs27qVkMeXY1a3YfrvPX0QJXSik3CPL34S/DOrHi94Po0iyU1pHBdf4aPnX+FZVSSv2sU2wo86f2dcvX1i1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyqXqd1FhE8oH95/n0SKDuz0WtH5rdGnbNbtfcoNndpZUx5r9mVK7XAr8QIpJyplmZ7UCzW8Ou2e2aGzR7fdMhFKWUsiktcKWUsik7FfgcqwNcAM1uDbtmt2tu0Oz1yjZj4Eoppf6TnbbAlVJKnUYLXCmlbMojClxE4kTkSxHZISLbRGS6a324iKwUkQzXMuy05zwgIj+KyC4Rucq69D/ncYrIRhFZ4bpvi+wi0kREFovITte/f38bZf+D6/tlq4i8LSIBnppdRF4TkUMisvW0dbXOKiK9RGSL67FZIiIWZf8/1/fMZhF5X0SaeFr2M+U+7bE/iYgRkUhPy10rxhjLP4BYINF1OwRIBzoDTwMzXOtnAE+5bncG0gB/4CJgN+C0+D3cC7wFrHDdt0V2IBm4zXXbD2hih+xAc2Av0Mh1fxEwxVOzA4OBRGDraetqnRVYB/QHBPgYuMai7FcCPq7bT3li9jPldq2PAz6l+qTCSE/LXZsPj9gCN8bkGmM2uG4XADuo/gEdRXXB4Fpe57o9ClhojCk1xuwFfgT61G/qfxORFsC1wNzTVnt8dhEJpfqb/FUAY0yZMeY4Nsju4gM0EhEfIBDIwUOzG2O+Bo7+YnWtsopILBBqjFlrqpvljdOeU6/ZjTGfGWMqXHe/B1p4Wvaz/JsDzATuA04/gsNjcteGRxT46UQkHugJ/ADEGGNyobrkgWjXpzUHMk97WpZrnVX+SfU3RNVp6+yQvTWQD7zuGv6ZKyJB2CC7MSYb+AdwAMgFThhjPsMG2U9T26zNXbd/ud5qt1K9ZQoenl1ERgLZxpi0Xzzk0bnPxqMKXESCgSXAPcaYk7/2qWdYZ8nxkCIyHDhkjEmt6VPOsM6qYzl9qP4T8yVjTE+giOo/5c/GY7K7xotHUf3nbjMgSERu+rWnnGGdpx5De7asHvceRORBoAJY8NOqM3yaR2QXkUDgQeB/z/TwGdZ5RO5f4zEFLiK+VJf3AmPMe67Vea4/YXAtD7nWZ1E9jvWTFlT/+WyFgcBIEdkHLAQuF5E3sUf2LCDLGPOD6/5iqgvdDtmvAPYaY/KNMeXAe8AA7JH9J7XNmsW/hypOX28JEZkMDAcmuoYXwLOzt6H6F36a6+e1BbBBRJri2bnPyiMK3LVX91VghzHm2dMeWgZMdt2eDCw9bf04EfEXkYuAdlTvaKh3xpgHjDEtjDHxwDjgC2PMTdgj+0EgU0Q6uFYNAbZjg+xUD530E5FA1/fPEKr3ndgh+09qldU1zFIgIv1c7/nm055Tr0TkauB+YKQxpvi0hzw2uzFmizEm2hgT7/p5zaL64ImDnpz7V1m9F9X1i3sQ1X+WbAY2uT6GARHAKiDDtQw/7TkPUr2neBceslcYuJR/H4Vii+xAApDi+rf/AAizUfZHgJ3AVmA+1UcQeGR24G2qx+rLqS6OqeeTFUhyvd/dwAu4zqa2IPuPVI8Z//Tz+rKnZT9T7l88vg/XUSielLs2H3oqvVJK2ZRHDKEopZSqPS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyKS1wpZSyqf8PUJabcVI7E1IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = 300.                  # Temperature, K\n",
    "Eps = .5                  # Emissivity\n",
    "ErrLevelEst    = 0.1      # Estimate of noise level, RU\n",
    "NoiseLevelTrue = 0.1      # RU\n",
    "bias = 0.\n",
    "noise = NoiseLevelTrue*np.random.randn(Nobs)\n",
    "yobs  = (1e3*plancknu(nu,T) + noise + bias) * Eps\n",
    "\n",
    "plt.plot(nu,yobs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"georgia\">\n",
    "Based on this \"observed\" spectrum, use inversion to get the best estimate for the cloud temperature and emissivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[294.88742172   0.46628455]\n",
      "[300.53535715   0.49983248]\n",
      "[299.91889972   0.50068808]\n",
      "[299.9155703    0.50070485]\n",
      "[299.91556656   0.50070487]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12364dfd0>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD8CAYAAABuHP8oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1hX5f/H8efNhykiQxBRlOFCHIDirGxalpY2LdPUTNv7197rm20rzb1ypZlpaTlSc4sCokxBQIYiQxFQ9of794d8+zY0GZ8P8/24rq7DWfd563X14njOue9baa0RQgjR+FjUdwFCCCFqRgJcCCEaKQlwIYRopCTAhRCikZIAF0KIRkoCXAghGqkqBbhS6lmlVLRSKkoptUIpZauUclFKbVFKJVQunc1drBBCiP+5bIArpdoDTwHBWuuegAG4F3gZ2Kq17gJsrVwXQghRR6r6CMUSsFNKWQItgJPASGBx5f7FwCjTlyeEEOJSLC93gNb6hFLqUyAVKAI2a603K6XctdYZlcdkKKXaXK4tV1dX7e3tXduahRCiWQkLC8vRWrv9fftlA7zy2fZIwAc4C3yvlBpb1QsrpaYAUwA6duxIaGholYsWQggBSqmUi22vyiOUG4BkrXW21roMWAMMBjKVUh6VjXsAWRc7WWs9R2sdrLUOdnP7xy8QIYQQNVSVAE8FBiqlWiilFHA9EAv8BIyvPGY8sM48JQohhLiYqjwDD1FKrQbCgXLgEDAHaAmsUkpN4kLI323OQoUQQvzVZQMcQGv9FvDW3zaXcOFuXAghRD2QnphCCNFISYALIUQjJQEuhBCNVJWegYtmpvQ8FJyCc5lQkAEFmRe2+w0HZ6/6rU0I8QcJ8OakpOBCGBdkVIbzqb/9fCG0y8/nk3W2FWfzHCjKt6PivBUo8On1Hq49A6H33eA/Clq41PefSIhmTQK8qSs5B1vfpXTnUk5nWZCXb09pgS36nAFDocKqsAKb4jJsi4uxLSnEtrgA25JztCOfdn9rqmy7DYkdM8ntMYuegS9j22Mo9L4HutwEVrb18scTojmTAG/Kjv1G5uznKV1RTIcTGXgAHn/arVEU27aiyNaBEht7cp09KGnhS4mDPeVOraho44xlBw9sO/tQmnYS25Ub8D52mE7J5zm/1ZmoTklUBLxBj+5PYuh5G/QeDV5XgIW8WhGiLkiAN0WFZzCuf5nI6WH4709A6Qoi+t5IkVdbtLsrVl4e2HfviltgAC7tPbEzGKrW7puvUpRzmvCX38B1y266Rx/EEGUkx8WL1C6HcQn+BW9fW+h114Uwd+9h3j+nEM2c0lrX2cWCg4O1DGZlZjHrSJ7xJjY/FNEuM4GTbbtz+sOX6DVh/OXPrabMsFBOvj4Vz4NhuJ0+ToWyIL2dP6f8HOnc5xgu3r4XHrH0uhsc25v8+kI0F0qpMK118D+2S4A3EQWnKF37LNEzjtMz9BDawoIjN9xG7zWLsLY1//PpY0uXUfLlArxiDtOy8DQl1vYc9+pBfk9NgH8C1t2ugVs+hdadzF6LEE2NBHhTpTVELCfy689psz4P9+xE0jx7UjZ7Kr63DK/zcoylpcS8+wF23/+CV1IUVuXF5Ld0J6lLF7xuTsX5rg8g8H5Qqs5rE6KxkgBvinJTKFj5BMfmnqXX4RDKLW2Jvv1OApfMxWBZ/683zp08ScKLb9Bm+17an4zjtHNHTg5vRa9bg2DENPkMUYgqulSAy+cCjVGFEfbP4sDj91L8ehRB4btJ7+hP3s7f6LtiYYMIb4CW7doRtHQ+7U/EcvjZF7AryqP7injCvkiifPoVkLSjvksUolGTAG9sso+SNe0WDj+8gr4rD2J/PpfDjzyJd1IE7v3713d1lxTw+ccU7NhMhkcX+u7fQ8ZXFqR+NgE2vwHlJfVdnhCNUsO4VROXZyzDuPNz9k3fiP+2RALOppHUtS8uPy0loJtffVdXJe79+2NMjuDQyDH02rSOsgV2HEzeTL+k3+HOeeDWrb5LFKJRkTvwxiAngeT/DCPmiY0M/nEnNqXniHvjTXyPhuLUSML7vwyWlgRtWMWx+XMptGtFv18PEPdFMbmf3wAH5114KSuEqBJ5idnQFZ5h75R78f85Gqf8kxzrPRDPjaux9Wj831UX5Zwm6ca76HHod3Id23HiFmd6jugBt02HljJ/qhD/JS8xG6MKI1HvP0jvH8OxKivk+LRpdD68r0mEN4Cda2t6hG8n/P9ewbqsmO4rj3Loy1TKvx4MCVvquzwhGjwJ8AYsa+VruM6PxaakgDMLZ+H99NP1XZJZ9PnkPxTs3Exahx4EHdhP1lcWpH4xBX55EcqK6rs8IRosCfAGqjxyLblv7aRtVjwJjz5Kh9Gj67sks2rbty8djoUSett9uOWk0mZ+DuFfboM518KpqPouT4gGSQK8IcqOJ/zh6XRL2Efs4Ovx/3pafVdUJwyWlgSvW86xhfPId3Cjz+YwEqaVkvflUNg/U15wCvE3EuANTXE+ex99nL4hv5Pu6UfX7b/Ud0V1rvvY+3GIDyUy+AY6JR6GGZroadNh06sS4kL8yWUDXCnVTSkV8af/8pVSzyilXJRSW5RSCZVL57oouEmrqCD6nYkEbThIQUtXHHdvxGBtXd9V1Qs7Fxd6HdxC+MuvYVFhxG/VUcK/2Azrn4WKivouT4gG4bIBrrU+qrUO1FoHAn2BQuBH4GVgq9a6C7C1cl3UQvYPb9FmbjgGYyn5i77BwUvmnwz+z7uc27WJHFcfArdEEPblHlj7KBjL67s0IepddR+hXA8kaq1TgJHA4srti4FRpiysuSmP+ZX8F3/B9XQK8U89Scfb76jvkhoMj6A+GHb9TGabzvTZGkbYV4fghwehvLS+SxOiXlU3wO8FVlT+7K61zgCoXLYxZWHNyplkjoyfSqfj4URfPYyen35U3xU1OK7d/LAN2cRJDz/67jhI2FfxsHIslBXXd2lC1JsqB7hSyhq4Dfi+OhdQSk1RSoUqpUKzs7OrW1/TV3qekEkPERi2m+Peven+20/1XVGD5ezljUPYb6R69qTvnhDCv06D5fdA6fn6Lk2IelGdO/CbgXCtdWbleqZSygOgcpl1sZO01nO01sFa62A3N+ke/RdaE/PmeAJ/CeGsUzta797QYIaCbahaebTDLWIHKV696ROyn0PTc2DJHVCcV9+lCVHnqhPg9/G/xycAPwH/nWhxPLDOVEU1Fzmr38dj9n60UhQsmYlDe8/6LqlRsGvtQtsju0jy7UNQ6D4OTy/AuHAEFJ6p79KEqFNVCnClVAtgKLDmT5unAkOVUgmV+6aavrymq/zob5x//nsc8zI4+uzTeA0fUd8lNSo2rVrhGbmLY137ERCxj5iZRozzh8G5i/5DUIgmqUoBrrUu1Fq31lrn/Wnbaa319VrrLpVLuf2pqrx0Yu59E6+0SI5cP4KAD9+v74oaJesWLfCJ2kO8/0B6Re4jbqYF5fNugrwT9V2aEHVCemLWtbJiQseNo3fEPo51DqbXrz/Ud0UXVVGhKSk3UlhaTn5xGcVlxvou6aIMVlZ0OryLuICr6BEbwrGZ1pTMGQa5x+u7NCHMTt6Y1SWtiXv9AXpvDCG7tTfuu34y+0vLvKIy5u9KYtvRLMrKNeUVFZRXaMqNGmPF/9aNRn1he+X633usWxssuLGHO6P7deCKTq5YWDScWeUNlpZ0Cd1GzKCh+If+TvzsfnQ03oLtQ+vAtUt9lyeE2UiA16GcHz7CY+ZOyi1tOL98Fm5tPcx2rYLiMhbuOc7cXUkUFJcz0NcFJydrDAaFpYXC0sICSwv113WDwmChsLJQGP60bmmhSDtTyNqIk6w/kkF7Jzvu6uvJ3cGeeDq3MNufoToMlpb4H9hG1JU30XPvFhLn9KVN2QgcJv8AbXvWd3lCmIXMyFNHyo/t5NSQKXhkJnLk9TcJeucNs1znfEk5i/ZeCO6zhWXc6O/OMzd0xb9dq1q3XVxmZEtMJqtC09h9LAeAKzu7cndwB270d8fWylDra5hC5HUj6LV9A8leQbiMz8dxyipo36e+yxKixi41I48EeF0oyCTmihH4R4YSevPdBP+yyuSXKCo1smT/cWbtSOLM+VKu82vDszd0pZeno8mvBZCeW8jqsHS+D03nxNkiHO2suD2oPfcEdzDJL4vaOnLLHfT+9UdSPHvTcvw5Wk9eCl6D6rssIWpEArweHbrzRoLWbCHObzBdonZiMJjuTrW4zMiykFRm/p5IzrkSruriynNDuxLUsW4Gh6yo0OxJzGHlwTQ2R2dSaqygV3tH7gn25LbA9jjaWdVJHRdz5Pb76L32O9La9cR63Hncn1gFnv/4f0CIBk8CvJ4URf9O+YC7KbGxxybuAA5uphkypqTcyMqDaczYfozM/BIGd2rNs0O70s/bxSTt18TZwlLWHjrBytB0YjPysbG04OaebbmnXwcG+rSulxefR8aMp9eKJZxs64flpArcn/lZXmyKRkcCvJ5EXDWEwN27OPjEM/T7+otat1daXsHqsHSmb0vgZF4x/bydeW5oNwZ1am2Cak1Da03UiXxWhaaxNuIEBcXl+LrZM/WO3vT3qftfMJEPPkzPhXNJ6RiAyyMVtHp8I7Qy3wtkIUxNArwe5O74DtubpnDGxZO2aZG1enRSbqxgTfgJvtqWQHpuEUEdnXh+aDeu6NwapRrOJ31/V1xmZGPUKT7fEk9abiFThvjy3NCu2FjW7QvPyJH30Oun74ntOphOT2isJ/8KtuZ5PyCEqV0qwKUjj7lozYlnvsSupIDMJ8fXKrwPJJ/hhs938OIPR3Cxt2bhxH6seXQwV3ZxbdDhDWBrZWBUUHt+ffoq7u3Xkdk7khg5fQ8xJ/PrtI5e61YRO/BausfvJWqhLcbl90F5SZ3WIISpyR24mZz67gtaj3uFtA498E0Kq3E7h1JzuX9eCG0cbHhtuD83dG/T4EP732yPy+LFH45wtrCUZ4d25eEhnTDU1bPxigqSuvfDNz6cfVdez6An28NdC8FC7mNEwyZ34HXJWE7eO8uwMJZT/OZTNW4m7lQ+ExYexLWlDaseHsRQf/dGHd4A1/q1YdMzQxjq787HG49yz+x9pJyuo/G8LSzoGL6Tk+26MmDPdnYvSIONL8tEyaLRkgA3g+Tpr9Hl6CHi/QbgP2H85U+4iOM55xk77wC2VhYse2gAbVrZmrjK+uNib82MMX348t5AEjILuPnLXSwPSaUu/jVoaW+PS9h2cp3bMWDrPvYu2Ae7a/9yWYj6IAFuamVF8MWvGA3W2H75do2aOHm2iPvnhWCsqGDppAF0cGkY3dVNSSnFyMD2bHp2CH06OvPqj5FMXHSQrHzzT5Fm27Ydhh2/UGzbkqANkRycuxwOLTP7dYUwNQlwE4t760l8UiKJCbgCn6FDq33+6XMljJ0fQn5RGd8+OIAu7g5mqLLh8HC049sH+/PObT3Yn3SaG6ftZMORDLNf16lnL86tXoLSFXRde4Ko2R9BwhazX1cIU5IAN6Wis7RauJMi21a4z/+s2qfnFZXxwIIDnMgtYv6EfmbrBt/QWFgoxg/2ZsNTV+HV2p7Hl4fzzHeHyCssM+t1PW4aRsaML2hReIY2qwtJnvU0pDePl+yiaZAAN6HIJybR7lQC0QOG0DYgoFrnFpaWM2nRQeIzC5g1rm+9dHipb53cWvLDI4N4bmhX1h/J4KZpO9mVYN6JsH0mTSL51ZdxPZOM5XeWZM0ZBznHzHpNIUxFAtxEjGfS8Vh7kPyWbei09JtqnVtSbuThJWGEp+YybXQQ13YzTXf7xsjSYMFT13dhzWODsbcxMG7+Ad5aF0VRqfkmlOj67jvEjZ9Ih5OR5C1zomD+nVBwymzXE8JUJMBNJHLSQ7ieSSP2+utw9uxQ5fPKjRU8810EuxJymHpHb4b3li7eAL09ndjw1FU8eIUPi/elMPzrXWb93NB/4Tziho6gS9J+kpa6Urr4TpnpXjR4EuAmUJp6GN/fDnHa2YuApbOrfF5FheblNZH8GnWK14d3555+VQ/+5sDWysCbt/qz/KEB5J4v5c6Ze4lMN1+o+m36icTAQQREbePQklbSW1M0eBLgJhA7/mlancsi4fabsW1ZtbGwtda8tyGG1WHpPH19Fx66ytfMVTZegzu7svrRwdhYGhg9Zx874s30XFwpOh3cSZp3D/of2MTepRp+fBgqKsxzPSFqqUoBrpRyUkqtVkrFKaVilVKDlFIuSqktSqmEymXdDEDdwBSEb6Hr/kOcbOtHv9lfV/m8L35LYOGe4zx4hQ/P3CDDm15OJ7eW/PjYYLxa2zNp0UHWhKeb50KWlnhE7CG7jReDdm5lx7cnYdMr0ltTNEhVvQP/EtiotfYDAoBY4GVgq9a6C7C1cr3ZOf7wm9gV55P+4F1VnqB43q4kvtqawN19PXl9ePdG3z2+rrRpZcvKhwfS38eF51YdZvaORLP03rR0dMThwO+cs3dhwOZ97F68H/ZMM/l1hKitywa4UqoVMASYD6C1LtVanwVGAosrD1sMjDJXkQ1VzsZv6RYRQbJXEP0/eK9K56w8mMr7G2K5pVdbpt7Zu0HN7t4YtLK1YuHEftwa0I4Pf43j3fUxVFSYPsTtvLxQv/1MhcFA4PoYQhasgMjVJr+OELVRlTtwXyAbWKiUOqSUmqeUsgfctdYZAJXL5vXtm9bkPP8VlsZS8l+YVKVT1h85yctrIrm6qxvTRgfV3Sh8TYyNpYEvRwfy4BU+LNxznCe/O0RJuek/M3Ts14/8FYuwLj1P57VZxMx9F9IOmvw6QtRUVQLcEugDzNRaBwHnqcbjEqXUFKVUqFIqNDvbvJ0y6lL6oo/pEneY+C4DCHj88csevz0ui2e+iyDYy5lZY/tibSnvj2vDwkLxxojuvHqLHxuOZDBhwUHyi03fc7PtyFFkfPoRTmfTcVxdwam5E+FsqsmvI0RNVCVF0oF0rXVI5fpqLgR6plLKA6BymXWxk7XWc7TWwVrrYDc3N1PUXP8qjJS+t4wKZUBNfeGyh4el5PLI0jD8PByYP6EfdtZ1OxtNU6WUYsqQTkwbHcjB42e4Z9Y+Ms0wGJbXU0+SMOVh2mdEk7+8FcWLR0NJgcmvI0R1XTbAtdangDSlVLfKTdcDMcBPwH/HSh0PrDNLhQ1Q4kf/h29yJLE9B9Pt9tv/9djiMiP/9/1h2rSyYfHE/rSyrb9Z2puqUUHtWTChH2lnCrnjm70kZp8z+TX8Zs0g9rpb6Jq4j9hFLTCuehAqzNc7VIiqqOq/458ElimljgCBwH+AqcBQpVQCMLRyvekrL8Fm+i8U2zjgPOuDyx4+e0cSyTnn+WBUL1q3tKmDApunIV3dWPnwIErKjdw5cy9hKbkmv0b339aT5NeXoIhtHJh/Gra8afJrCFEdVQpwrXVE5WOQ3lrrUVrrXK31aa319VrrLpXLM+YutiGIefYhPE/GE933SjoMHPSvxx7POc+M348xvLcHQ7o2kcdHDVjP9o788OhgnOysuH/efrbGZpr2AkrRMXQHGW270H/XbkJmbofQhaa9hhDVIG/SqqMoD7cVOymwd6Xjt1/966Faa978KRprgwVvjvCvowKFV2t7Vj86mK7uDkz+NpTvDpj2haOlvT0tdv9CgYMrARsjOTJ9GiT9btJrCFFVEuDVcOTBcbidTiXmqmtw69T5X4/9NeoUO+OzeW5oV9yb0HRojYFrSxtWTB7IVV3ceHlNJF9tTTBphx/HTp0p+G4BFcqA9w+nSJ3+GOQkmKx9IapKAryKjDlpeG0I4YxTB/yXzfrXY8+VlPPuzzH4e7TigUFedVSh+DN7G0vmjQ/mjj7t+XxLPK+vjTJph58Ow4aR/MHb2BXlYbukhLw5o6GwWTxFFA2IBHgVRY+dgGNBFkdvvgEHl9b/euwXW+LJLCjmg9t7YmmQv+L6YmWw4LO7A3j0mk4sC0nllTWRJg3xHs89w5GJD+GWk0zu7DLKl42B8lKTtS/E5Ui6VEHxsXA67wznVJuu9F3w75M1xJzMZ9He49zXvyNBHZvl+F4NilKKl4b58dR1nVkZmsYb66JM+jil79zpHLpuJN6pUcR/mQMbnpWBr0SdkQCvgoQJT9Oi6CzJ996Kte2ln2dXVGheXxuJk50VL93kV4cVist5dmjXP+7E3/4p2qQhHrDpe2L9r8Q/OoyIT/bD3qqPSilEbVRt+LxmzJiXhU9ELCc8/On/+Uf/euzK0DTCU8/y6d0BOLaQDjsNiVKKF2/qRrmxgrm7krE0WJhsJEiDpSXee9aT2msIvfdGEDW1hJ6fdAK/4SaoXIhLkzvwy4h5+nFanj9N2lV9MBgu3QX+9LkSpv4aR38fF+7s074OKxRVpZTi1Vu6M2GwN/N3JzN1Y5zJ7sTtnByx2bSc0y5edNuYSNLU5yDjsEnaFuJSJMD/jda4bD7MOfvW+E//7F8P/fDXOM6XlPP+qJ4yvncDppTirVv9GTuwI7N3JPHZ5niThbi7fw9y5n5GsY0D7quyOP3FGMjPMEnbQlyMBPi/SFv8Oe0zEjjmH0grt0uPlnsg+Qyrw9J56Cpfuro71GGFoiaUUrx7W0/u7deB6duP8eVW033D3f3224l59QUsy0tRc3MpmXc3lBaarH0h/kwC/F+c+3wlRmXA7tXHLnlMmbGC19dG0t7Jjqeu//fOPaLhsLBQ/Of2XtzV15NpvyUwY/sxk7U94NUXOHDfRBzzszj96UlYPUXm1RRmIQF+CcXJ0fgejSHFK4Buo+645HHzdycTn3mOt2/rQQtreSfcmFhYKD66sze3B7Xnk01Hmb0j0WRtX7V4OiHX3E27jGSS3jkA2y8/8JkQ1SUBfglxTz6PTel5sm8dcslj0nML+fK3BG7o7s5Qf/c6rE6YisFC8cldvf+Yom3eriSTtR38yyIiAobhG59A/GtL4PBKk7UtBEiAX1x5KZ77Yjjj1IGgjz+85GHv/BwDwNu3yWBVjZmlwYIv7gng5p5teX9DLIv3HjdJu9Y2NnhtWsgxn0F0DUkg6Y3XIO2ASdoWAiTAL+ro+y/ieiaNxKCAS3bc+S0mky0xmTx1fRc8nVvUcYXC1CwNFnx1XxBD/d1566doloWkmKRdZ/e2sOILTrl1o8Ov6WR+MA7OppmkbSEkwC/CcunvlFq1wPOzty+6v6jUyNs/R9OlTUsmXelTt8UJs7EyWDB9TBDX+bXhtR+jWHnQNEPRdh4wgLQPX6HExgH7pTkUTr8LSkw/a5BofiTA/yZ37y94J0eT2CkQj6C+Fz3m620JpOcW8f6onjI5cRNjY2ngm/v7MKTrhaFoV4elm6TdfpPGEzZ5MrbF5yj4OhPjionyZYqoNUmfv0l78SMMFeUUPTjqovuPZRUwd1cSd/bxZIDvv49KKBonWysDc8b15YpOrryw+jDrIk6YpN2rv/gP+2+8B/esFJI/jIKt75ikXdF8SYD/ibHgDD6HYzjp7kfAc8/9Y7/WmtfXRtHC2pJXb5HBqpoyWysDcx8IZoCPC8+ujGD9kZMmaXfQz4s43PtGOsfHEfXuGohYYZJ2RfMkAf4n0c8+gcO5HNKuCLrouCc/HjrB/qQzvDism0xQ3AzYWRuYP74ffb2cefq7CH6Lqf0cmwaDAZ9ty0jxDMR/TwxxH7wHqSEmqFY0RxLg/6U1zhsPcb6FC37TP/3H7rzCMj7YEEtgByfu69exHgoU9cHexpKFE/vTs10rHlsezr7E07Vus1VrV1g9k7OO7fH+KZWTn0yEs6adu1M0DxLgldKXz8DzxFGO+QXg6NHuH/s/3hRHbmEp74/qiYWFDFbVnLSsDHEvlxZM/jaUI+lna92m14CBpH38FigL7JbkkT/rHigpMEG1ojmpUoArpY4rpSKVUhFKqdDKbS5KqS1KqYTKZaOefib/k6VopbB95eF/7ItIO8vyA6mMH+xNz/aO9VCdqG8u9tYsmTQApxZWjF9wgGNZtQ/bgMmTiJg0Gcf8U5ydnofxuwehwmiCakVzUZ078Gu11oFa6+DK9ZeBrVrrLsDWyvVGqTg1Ht+4aFI69KbbXaP/sf/TTUdpbW/Dc0O71kN1oqFo62jL0kkDMFhYMHbeAdJzaz/K4MAZnxF23Sg6nogj9vNk+TJFVEttHqGMBBZX/rwYuPh3d41A7BPPYltyjqzhV/5j35H0s+w+lsPkq3xwsJVZdpo7b1d7lkzqT2FpOWPnhZBdUFLrNvtsXMlRv8H0PBJG2Ee/wqFlJqhUNAdVDXANbFZKhSmlplRuc9daZwBULi89YHZDVmHEc18MuY7tCfr8k3/snrUjEQdbS8YMkBeX4oLuHq1YOLE/mfklPLDgAHlFZbVqz2BpSfudaznVpgsBvx8h8pOPIWWfiaoVTVlVA/wKrXUf4GbgcaXUpYfo+xul1BSlVKhSKjQ7O7tGRZrT0f+8glvOcRID/znuSVL2OX6NOsUDg7zk7lv8RV8vZ+Y80JdjWQU8uOgghaXltWqvpZsbxnWLKbJzxGddOinTHoZc04zHIpquKgW41vpk5TIL+BHoD2QqpTwAKpdZlzh3jtY6WGsd7ObmZpqqTcjw7W+UWdri8elb/9g3Z2cS1gYLJgyW8U7EP13VxY0v7w3iUGoujy4Np7S8dl3j2w8cRNrH72BTeh67pefJm3sfFOebqFrRFF02wJVS9koph//+DNwIRAE/AeMrDxsPrDNXkeaSG7IF76RoknwDaR/c/y/7TuUV80N4OvcEd8DNQTrtiIu7pZcHH97Rix3x2Ty7MgJjRe3m1/R/9FGiJk6iTU4SWbNLKV85Ub5MEZdUlTtwd2C3UuowcADYoLXeCEwFhiqlEoChleuNSuqLH2JpLOX8hBH/2LdgTzIVGqYM8a2HykRjMrpfR167pTsbIjN4fW1krSdJDpozg8grb6JLUhiHZ+TAljdNVKloai47B5jWOgkIuMj208D15iiqLhgL8/E+FM2pNl0JePGvX0DmFZaxbH8KI3p70MFFxvoWlzd5iC95RWVM336MVnZWvHJz91q157/tZ45360dQ6G72fqkZ7LYE+owzUbWiqWi2PTGjn3sSx4IsUgYF/GPck2/3Hed8qZCFIMUAACAASURBVJFHru5UP8WJRun5G7sybqAXs3ck8c3vtZsk2WBlhfveTeQ6taPv1nBCp30FybtMVKloKpptgDttCKXQ1gm/GZ//ZXtRqZGFe49zbTc3unu0qqfqRGOklOKd23owMrAdH288WutZfezaumNcu4wKCwOd150k4esnIDveRNWKpqBZBviJVXPwPBF3YdyT9p5/2bcqNI0z50t59JrO9VSdaMwsLBSf3h3A9X5teH1tFD8drt0wtG2GDCH9w3dxOJdFy+WlZM2+D87nmKha0dg1ywDP/XARADYvTf7L9jJjBXN2JhHs5Ux/H5d6qEw0BVYGC2bc34d+3i48tzKC7XEX/cK2yro88zRxkx/GIzOO/AUGir+9F8qKTFStaMyaXYCXZCTTKS6aNM+edLv3/r/sW3/kJCfOFvHoNfLsW9SOrZWB+eOD8fNw4JGlYRxIPlOr9nrMmk700FvpnHyQozPLMP7wsEzJJppfgMc8/jR2xflkDhv4l+0VFZqZvyfSzd2Ba7s1zlEBRMPiYGvF4on98XS2Y9Kig0SdyKtVez02reNYjwEERO4kdHo6bHvXRJWKxqp5BXhFBe12R3O2lQeBX077y67tR7OIzzzHI9f4ynjfwmRat7RhyaQBtLK7MAxtYnYtZqNXCp/Q3znR3o9+e3ew94utELb48ueJJqtZBfjRT9/APTuJxN69sbaz+8u+mb8n0t7JjhG9/zmZgxC10c7JjqUPDUApGDcvhBNna/782mBri/OB38h1akfwljBCP50BidtMWK1oTJpVgKv5Gykz2OD+yRt/2X7w+BlCU3KZMsQXK0Oz+isRdcTH1Z5vHxxAQUk54+aFkHOu5sPQtmjXHuMvqyi1ssPvpyRiP3sesmJNWK1oLJpNWp09tAufpGiSfQPwHHjFX/bN/D0RF3tr7gnuUE/ViebAv10rFk7ox8m8Ih6YX7thaNsMHETmrGlYlxXS5rsznPxqLBTUftJl0bg0mwA//tw7WJWXcO7+m/+yPTYjn21xWUwc7I2d9T9nohfClIK9XZg9LpiErAIeWnyQotKaD1TVadw44l96AaezJ6hYUEzBgtFQWvtZgkTj0SwC3FhciHd4DFmuvgS8/tfHJ7N2JGJvbeCBQd71U5xodq7u6sa00UGEpeTyyNKwWg1D2/O9d4i8+348M2LImH6e8lUyr2Zz0iwCPPr/nsQpP4OUgUF/Gfck9XQhPx8+yZgBHXFsIRM2iLozvLfphqENXLmYqEFD6XoslKjP02X0wmakWQS4w4ZQim0c6Dz9s79sn7srCYOFYtKVMmSsqHt/Hob2tR9rNwxt952/kNi5L4GH9hDy0Q44MNeElYqGqskHeO7hPXRIi+O4bw+cvbz+2J5dUMKq0DTu7ONJW0fbf2lBCPOZPMSXJ67tzHcH05j6a1yNQ9xgaYln6FYy2nYleEcIoR/OhfjNJq5WNDRNPsCTX/kAS2MpRXdc+5fti/YmU2qskAkbRL17/sauPDDIi9k7k/jm98Qat2Pj6Ij11jXkO7Sh1y8xRE19CU5FmrBS0dA07QDXGo/wY+S18qD32//rdlxQXMa3+1K4uWdbfN1a1mOBQlwYhvbtW3swKrAdn2w6ypL9NR+GtrV/D/KXz8VosMLr+xOkfjoe8ms3IqJouJp0gKesmo9HZgLHu/lhsPzf5EPLQ1IpKC6XCRtEg2Fhofjk7gBu6N6GN9dFsS7iRI3b8rrlFpI/eBvb4nzsFuWSO3M0lNSiC79osJp0gJ/9YikAtk+N/2NbcZmRebuTubKzK709neqrNCH+wcpgwfQxfRjg48Jzqw6zNbbmHXN6PPs0kZMfofWZVAqm51K6fAIYy01XrGgQmmyAG0uK6RiXwCm3znQb+78A//HQCbILSmTIWNEg2VoZmDe+Hz3ateKxZeHsTzpd47b6fDONiJvuoGN6NEkfpWBc+yTUcsJl0bA02QCPnfoGznknORHk98c2Y4Vm9o5Eens6MrhT63qsTohLa2ljyaKJ/eno0oKHFodyJP1sjdvq8+v3RAdejd/RUKI+PARbZQjapqTJBrjF97swWljh+eH/el7+GpXB8dOFPHp1J5SSIWNFw+Vib82SSQNwanFhGNqEzIIat9Vt/yYSOwUTEBZC+Du/wL4ZJqxU1KcqB7hSyqCUOqSUWl+57qKU2qKUSqhcOpuvzOopyTmJd2IcaZ7dce/THwCtL0zY4Otqz4092tZzhUJcXltHW5ZOGoDBwoJx8w+Qerpm45xY2tjQ7sBmUj17ErT7IGFvLYLD35m2WFEvqnMH/jTw5zErXwa2aq27AFsr1xuEqBdepEVxHjlXB/6xbVdCDtEn83n4al8MMmGDaCS8Xe1Z+lB/isuN3Dd3f43HErdzccZh989ktulE4NYIDr35McRvMnG1oq5VKcCVUp7AcGDenzaPBP47HchiYJRpS6s55x3RFNs40OOLz//YNvP3RNq2smVUUPt6rEyI6vNr24qlkwaQX1zGmLn7OZVXXKN2nL28UVtWc9axHT1/iSXqrZcgNcTE1Yq6VNU78GnAi8Cfh01z11pnAFQuG8REkjlhu+iQGkuKT3fsWl94UXkoNZd9Sad56CofbCxlyFjR+PRs78i3D/bn9LlSxszbT1ZBzULcvXcA51YvoNDOkc4/JhH/3mTIjDFxtaKuXDbAlVIjgCytdVhNLqCUmqKUClVKhWZnZ9ekiWo5/sZHWBlLKLln6B/b5u9OppWtJff272j26wthLkEdnVk4sR+n8ooZOy+E0zWc1cfruuvJmPslRgtL2i0/QerUMXA21cTVirpQlTvwK4DblFLHge+A65RSS4FMpZQHQOUy62Ina63naK2DtdbBbm5uJir7ErSm3aFj5Dm40+ONt4ELg1Ztij7FXX070NLG8t/PF6KB6+ftwrzxwaScLmTs/AOcLSytUTt+o0eTOPVdrMpLcFyQSdand8L5HBNXK8ztsgGutX5Fa+2ptfYG7gW2aa3HAj8B/+0hMx5YZ7YqqyhpxRw8TsWT4tf9j67zq0LTKDNq7h8od9+iaRjcyZW5DwSTmHWOcbWYmq33U08Q9cILtDifi8XMLHJn3AklNf9cUdS92nwHPhUYqpRKAIZWrtervK9XotC0fH4ScKHjzooDqQzybU0nGbRKNCFDuroxc2wf4k7lM2HhAc6V1KybfN/33+LQpIdxzj1B8ecnKVpwL5TXfMJlUbeqFeBa69+11iMqfz6ttb5ea92lcnnGPCVWjbGkGK+4BDLdfPEdPRaAnfHZpOcWMXag12XOFqLxub67O1/f14cj6XlMXHiAwtKahXj/2V8SeudY2mYmkv2fZMpXyLRsjUWT6YkZ+cFruJxN52Sf7n9sWxaSgmtLG4b6u9djZUKYz7CebZk2OpCwlFweWhxKcVnNgnfA94sIu24UHdNjOf5OFMZ1z8q4KY1Akwlw6zX7MFpY0mHqOwCcOFvEtrgsRvfzxNqyyfwxhfiHWwPa8endAexLOs2UJWE1DvGgzd9zuN+NdE48Qvwbu2D7hyauVJhak0i2olNpeCfFkt7eD9fAvgB8dyAVDdwnnw6KZuCOPp5MvaMXO+OzeXxZeI1mujcYDPTcvZ5Y/6voHh1B5Es/yNyaDVyTCPCol16mRdFZzlx3IbzLjBV8dzCNa7u1wdO5RT1XJ0TdGN2vI++N7MHWuCyeWnGIMmMNQtzaCt+QDST6BtMrNJzD/zcTon4wQ7XCFJpEgLvsiqXYuiXdK7vO/xaTSXZBCfcPkLtv0byMG+TNGyP82Rh9iudWHcZYUf3n2DYtHXDbt5609j3pvSuCwy+8D4nbzFCtqK1GH+BZB7bRMTWGFN/u2Dq7ALA0JIX2TnZc061B9O4Xok5NutKHl4b58fPhk7yw+jAVNQjxVm3csd32A5ltOtNzcxRRLz4HqfvNUK2ojUYf4Mff+gIrYwnl990MQFL2OfYcO819/TvIqIOi2Xr0mk48e0NX1oSf4NUfI2sU4m5du1K2bjFnHdvR7aejxL40WQa/amAad4BrTfvDieS3bIPfqxcmblhxIBVLC8U9/TrUc3FC1K+nru/M49d24ruDabz1UzS6Bp8Fdhg4iNPfzuB8Cxc6rz5G7EsPQdoBM1QraqJRB3j8tzNolxFHSvcLs84Xlxn5Piydm3q0pY2DbX2XJ0S9Ukrxfzd2Y8oQX5bsT+H1tVE1uhPvOmIEmd9+w3n7yhB/YRKkHTRDxaK6GnWAn5u5FoXG8cVHAPglMoOzhWXy8lKISkopXrnZj0eu7sSykFT+b/VhymvwdUq3kSPJXPwN5+xb0/mHY8S98CCkh5qhYlEdjTbAjaXFeMUnkOXqQ8e77gNgWUgqvq72DJIJi4X4g1KKl4Z147mhF56JP/1dRI2+E+82ciRZ314I8U4/JBD3/ERIr9Eo08JEGm2AR7zzKq1zUznZtwcAsRn5hKXkMmZAR5mwWIi/UUrx1PVdeO2W7myIzODRpTXrsdntttvIXDKTc/ZudFqTQNzzE+CEhHh9abQBbrvuAEZlwPvTD4AL457YWFpwV1/Peq5MiIZr8hDfPzr7PLQ4tEYDYPndeitZS2ZS0NKNTj/EE/fcBDgRbvpixWU1ygAvOJmEd3IsJzy74dSzN+dKyvkx/AQjerfDqYV1fZcnRIM2bpA3n9zVm72JOYxfcICC4uqPJ97t1hFkL5lFgYP7hRB/djycPGSGasW/aZQBHv3SG9gXniH3hn4ArIs4wflSo0zaIEQV3R3cgS/vDeJQ6lnGzgup0cw+3UYMJ3vJTAoc2tJ5TTxHn3lAQryONcoAd90bT4l1C7p/Pg2tNUv3p9LdoxVBHZzquzQhGo1bA9oxc2xfYjMKuHfOfnJqMMdmtxHDyV46k7xWbem0Jp6jTz8AJyPMUK24mEYX4Bl7NuGVGkmqb3esnZyISDtLbEY+YwfKy0shqmuovzvzJwRz/PR5Rs/eR2Z+9We77zb8FnKWVIb4j0c5+vQ4yDhshmrF3zW6AE957xusykuoGHcrAEv3p2JvbWBkYPt6rkyIxumqLm4sntifU3nF3DN7H+m5hdVuo9vwW8hZOou8Vu0uhPhTYyHjiBmqFX/WuAJcazyPJFLQ0pVuL7/B2cJS1h85yaig9jLjvBC1MMC3NcsmDyT3fCn3zNpHcs75arfR7ZabyVk2i7OO7StD/H4JcTNrVAEeteBr2p2KJdW/O1hY8EP4CUrKK7h/gMx5KURtBXZwYsWUgRSXV3DP7H3EZ1Z/hvpuNw/j9NKZnHVsT+c1cRx98n44FWmGagU0sgAvmbceC12B86tPobVmWUgKfTo64d+uVX2XJkST0KOdI6seHogCRs/eR9SJvGq30e3mYeQsm0WuUwc6/xjH0SfGwKko0xcrLh/gSilbpdQBpdRhpVS0Uuqdyu0uSqktSqmEyqWzOQstLy7E62g82a5etBt5F/uSTpOUfV7uvoUwsc5tHFj18CBaWFty39z9hKfmVrsNv2E3kbNsJmecL4R4/OP3ySeGZlCVO/AS4DqtdQAQCAxTSg0EXga2aq27AFsr180m/O3XcM1NIaNfTwCW7U/FqYUVw3t7mPOyQjRL3q72rHpkEK3trRk3L4T9Saer3YbfsJs4vXQmZ5w70mltHHFTxkHidjNU23xdNsD1BecqV60q/9PASGBx5fbFwCizVFjJbkNYZdf5qWQVFLMp+hR39fHE1spgzssK0Wy1d7Jj1cODaOdkx/gFB9gcfarabfgNu4nTy2eS7epLtw0xxDz4FESuNkO1zVOVnoErpQxKqQggC9iitQ4B3LXWGQCVS7PNX5aXmoBvUjQnPbvQyr8n34emU16hGSPDxgphVm1a2fLdlIF092jFw0vDWLQnudpt+N14I2W/LifNMwD/nVHETngH9s00Q7XNT5UCXGtt1FoHAp5Af6VUz6peQCk1RSkVqpQKzc7OrlGR0a+8jX3hGfJuHIixQrM8JJUrOrfG161ljdoTQlRd65Y2rJg8kKHd3Xn75xjeXx9T7YkhOvTpS4vdP5LQaRDdw6NImDgD4y+vQw1mCRL/U62vULTWZ4HfgWFAplLKA6BymXWJc+ZorYO11sFubm41KrLNgThKrezw+2waO+KzOHG2SF5eClGH7KwNzBzblwmDvZm3O5nHl4dXezhaVy9vPMM3EdnrerocjSb90R8oXTEFjNUfEVFcUJWvUNyUUk6VP9sBNwBxwE/A+MrDxgPrzFWkz46fyZ3xBZaOjizbn4qbgw1D/d3NdTkhxEUYLBRv39aDN0b4szH6FGPm7ufM+eoNgmXXyoHuYb9ycPBteKXGcfq5bRTMvhtKq9/7U1TtDtwD2K6UOgIc5MIz8PXAVGCoUioBGFq5bhaGdu1wn/ww6bmFbDuaxb39OmBlaFSfsAvRZEy60odvxvQh+mQ+d3yzh+PV7LVpaWVFvz3r2DdiDG2yUyh5I4ycj4dD4RkzVdx0VeUrlCNa6yCtdW+tdU+t9buV209rra/XWnepXJr9b/+7A2ko4N7+8vJSiPp0cy8Plk8eSF5RGbd/s4ewlOp/Kz7o52WEjH8Ex/wsrD6OIe31WyAv3QzVNl2N5ja2zFjBdwfTuM6vDe2d7Oq7HCGavb5ezqx57Aoc7awYM3c/v0ZmVLuNwQumc+j/XsKmtIjWM+OIf34kZMWaodqmqdEE+OboTHLOlcjLSyEaEB9Xe354dDD+7Vrx2PJw5u+u/meG/T98h7jPPsJosML721giHx8LqSFmqLbpaTQBviwkhfZOdgzpWrMvWYQQ5vHfzwyH9WjLe+tjePunaIzV/Mww8PFHOfntbM7bt8Z/zREOPfw4HP3VTBU3HY0iwBOzz7E38TRjBnTEYCGTNgjR0NhaGZgxpg8PXenDor3HeXRpGEWl1fvMsNsdd3D+11XktPYmcGM4hya9BeFLzFRx09AoAnx5SCpWBsU9wR3quxQhxCVYWCheH+HPW7f6syU2k/vmVn+aNs+Bg7AK2UiaZ0+C9hzkyKTpGLd9JB1+LqFRBPhNPdry8s3dcXOwqe9ShBCXMfEKH2aN7UvcqXzu+GYvSdnnLn/Sn7j4dMLt8A6OdelP74gDJDz2PaVrn4eK6t3RNweNIsD7+7gw6Uqf+i5DCFFFN/Voy4rJAzlfUs4dM/dy8Hj1vjK2c3HGO2oX0UHX4Hc0jBPPbCVv9p1QVP3PFZuyRhHgQojGJ6ijM2seG4xzC2vunxfCqtC0ap1vaW1Nj/DtRFx7Gz6pRyh94zCJr94snxn+iQS4EMJsvFrbs+bRwfTzdubF1Ud4Zc2Rao+hErhtHWGPPIljfibtZsUS9sgEiPnJPAU3MhLgQgizcra3ZvHE/jx6TSdWHEjjntn7SM+t3tgnfWd+xfEFcyixaUnQ2lDCH5yKcdPbzf65uAS4EMLsLA0WvDTMj9nj+pKcfZ5bv97NroTqDS/ddexYKkK3k9ahJ31C9nH8kXXkfnMXFJ01U9UNnwS4EKLO3NSjLT89eSVtHGx5YMEBZmw/Vq2xxV06d8XzWCiHr7wZ3+OHMb4VTsIrN0NWnBmrbrgkwIUQdcrH1Z4fHx/MbQHt+GTTUaYsCSWvqKzK5xusrAjY9QsRTzyHw7nTdJgVTejkCRC73nxFN1AS4EKIOtfC2pJpowN5+1Z/fj+azW3TdxObkV+tNoK+/pTUpQsptHOkz88HiXjgPxg3vwsVFWaquuGRABdC1AulFBOu8OG7KQMpKjVy+zd7+PFQ9YaT7XL33ajwnaR49SYwdB+pU37gzDd3QXGemapuWCTAhRD1KtjbhfVPXUlvTyeeXXmYt9ZFUVpe9btoZx8fOiYcJOKaEXilRsEbBzn64s2QfdSMVTcMEuBCiHrXxsGWZQ8NYPJVPizel8K9c/ZxKq+4yucbLC0J3P4zh597kRaFZ/GefYTQSRMg7hez1dwQSIALIRoEK4MFrw33Z8aYPsSdKmDE17vYl3i6Wm0EffohGauWcM6+NX1/DiFizHuUb36/yT4XlwAXQjQow3t78NMTF2b6GTs/hDk7E9HVGI3QZ+QorCL3kNSpD4GH9nPioVXkTL+nSY6jIgEuhGhwOrdxYN0TV3JTD3f+80scjy0Lr9anhq3ae+Idt5+IG0bRIT0aw5v7iHziZkjaYcaq654EuBCiQWppY8mMMX147ZbubI7JZNi0ndXqvWmwtCRwy48ceeV1bEoL8V8UTsTolyld9yKUVf35ekMmAS6EaLCUUkwe4suPjw2mhbWBcfMP8Oa6KApLy6vcRuAH71CwYzNpHfwJDN3PmclrOPbiDXAqyoyV143LBrhSqoNSartSKlYpFa2Uerpyu4tSaotSKqFy6Wz+coUQzVFvTyc2PHUVk6704dt9Kdzy5S7CUqr+TNu9Xz86HAsldPR4nPNO4TU9jPC7H8K444tG/YKzKnfg5cDzWuvuwEDgcaWUP/AysFVr3QXYWrkuhBBmYWtl4I0R/qyYPJAyo+buWXv5eGNclb8ZN1haEvzdIjJ+XEW2mzd9dodw8v65nHj/FjhbvbHKG4rLBrjWOkNrHV75cwEQC7QHRgKLKw9bDIwyV5FCCPFfgzq1ZuMzV3FXX0+++T2RkTP2VKsbvvctt9AmOYLwG27H41Qirf+zh0P33gVHvjdj1eZRrWfgSilvIAgIAdy11hlwIeSBNqYuTgghLsbB1oqP7wpg3gPBZBcUc9v03cz8PRFjFUc2tLSxoc+WNRydNYuClm4EbTpA4qj3yJ1xX6P63LDKAa6Uagn8ADyjta7yrzul1BSlVKhSKjQ7u3rj/wohxL+5wd+dTc8M4Ybu7ny0MY57Zu/jeM75Kp/f46GJtIwP5fCAm/A5Ho3Vi5s58sBNjeZzwyoFuFLKigvhvUxrvaZyc6ZSyqNyvweQdbFztdZztNbBWutgNzc3U9QshBB/aN3Shm/u78O00YHEZxZw85e7WLo/pcqdf+xcXAjYv5GI9/5DmZUtvVcfJG74MxR9/3yD/9ywKl+hKGA+EKu1/vxPu34Cxlf+PB5YZ/ryhBDi8pRSjApqz+ZnhxDs7czra6MYv/BgtcZT6fPaS+iIXcT0HIJfXAQlDy0n7pHrG/Tnhupyv6WUUlcCu4BI4L+ve1/lwnPwVUBHIBW4W2t95t/aCg4O1qGhobWtWQghLklrzdL9KfznlzisDIr3RvXktoB2XLgXrZqQJ1/Af+Ei7AvPEBsQQLdP78fy2mfAwmDGyi9NKRWmtQ7+x/bqjDFQWxLgQoi6kpxznudXRRCeepbhvTx46zZ/2jjYVvn89IMHKbrvSbokhpDj0oG80V50em0atO9rxqov7lIBLj0xhRBNko+rPd8/MpgXh3Vjc8wprvt0B/N2JVFmrNp34579+uF7dA/77n2IFoUF+M7cQ8J1E8ifOwkK//VhQ52RABdCNFkGC8Vj13Rm0zND6OvlzPsbYhn+1S72JuZU7XyDgUEr5pK17RdielxNp4QYbJ5YSdSw6zDunVfvvTglwIUQTZ6vW0sWTezHnHF9KSw1MmZuCE8sDycjr6hK53sPGkSPqO0c/PBjclp703PXYfKGv0XSlCsh47CZq780eQYuhGhWisuMzNqRyMzfE7FQiiev78ykK32wsazaC8qS8+cIHfckvTdvwOF8Nse9u+H04hU4TfgM7JzMUrO8xBRCiD9JO1PIe+tj2ByTiY+rPW/d6s813areoTwlLJTcya/T88h2tFIk9u1Cl8+fwzBoAlTji5eqkJeYQgjxJx1cWjDngWAWTewHwISFB5nybShpZwqrdL5X32ACwzcS+tnnnGjXC7+QKM4Pe5Hj4wdDZrQ5S/+D3IELIZq9knIj83cn8/XWY1RozaPXdOKRqztha1XFxypFReyb9BwBv6zHOS+dE+19aPniNThOmQa2rWpdnzxCEUKIy8jIK+KDDbGsP5KBp7Mdb47wZ6i/e5U7AaUeOcypKW8RGL4Fg7GEpMCu+H71EobBD9TqsYo8QhFCiMvwcLRj+pg+LJ88ADsrA1OWhDFx0UGSqzhAVsfeAfTfv5awr78k0WcQXcJjKR76NMfvHQjZ8SavV+7AhRDiIsqMFXy7L4VpW+IpLjdyX/+OPHZNZ9o6Vq03Z0lREbsffYVeGzbQJucYh8aOJWjJkhrVIo9QhBCiBrIKivliSwLfh6ZhYaEYO8CLR67xrXK3/LToaJKf/pDu89/Hzcu7RjVIgAshRC2kni7k620JrDl0AiuD4oFB3jw8xJfWLW3Mfm0JcCGEMIHknPN8tTWBdREnsLUyMGGwN5Ov8sXZ3tps15QAF0IIEzqWVcCXW4+x/shJ7K0tefBKHyZd6YOjnZXJryUBLoQQZnD0VAHTfovn16hTONhaMvkqXyZe4Y2DremCXAJcCCHMKPpkHtN+S2BLTCZOLayYMsSX8YO8sbexrHXbEuBCCFEHjqSf5Yst8Ww/mk1re2seuboTYwd6YWdd89l8JMCFEKIOhafm8sWWeHYl5ODa0oav7gtkcCfXGrV1qQCv/b29EEKIf+jT0ZklkwZwIPkM07cfw8fV3uTXkAAXQggz6u/jwrc+/c3StoyFIoQQjZQEuBBCNFKXDXCl1AKlVJZSKupP21yUUluUUgmVS2fzlimEEOLvqnIHvggY9rdtLwNbtdZdgK2V60IIIerQZQNca70TOPO3zSOBxZU/LwZGmbguIYQQl1HTZ+DuWusMgMrlJWcCVUpNUUqFKqVCs7Oza3g5IYQQf2f2l5ha6zla62CtdbCbm5u5LyeEEM1GTQM8UynlAVC5zDJdSUIIIaqiph15fgLGA1Mrl+uqclJYWFiOUiqlhtd0BXJqeG59k9rrx/+3d3ahVlRhGH5ePGVqiZVE5bG0iMCrtC60IiKjX9GrwEgsqvt+iFKEoksrKiIoQgkqU8RMQoi8qMsyzPKn1FK0PKblTSUFpfR2tOi+JQAABElJREFUsdap6bD30V3HM2vB98BwZr61B57ZzHx7Zq0136nVvVZvCPfTxaWdgiethSJpNXAj6eB+AJ4CNgBrgUuA74C7bA8d6BxRJG3pVAugBsK9HWp1r9Ubwn20OekduO27uzTNHWGXIAiCoAfiTcwgCIJKqSmBv9a2wP8g3NuhVvdavSHcR5VRrQceBEEQjBw13YEHQRAEDYpI4JKmSvpI0i5JX0p6KMe7Fs2StFTSXkl7JN3anv3fPmMkfS5pY96uwl3SJEnrJO3O3/+citwfyefLTkmrJZ1VqnuvReG6uUq6WtKO3PaSJLXk/mw+Z7ZLelfSpNLcO3k32h6TZEmTG7EivHvCdusLcBEwK6+fA3wNzACeAZbk+BJgeV6fAWwDxgLTgX3AmJaP4VHgbWBj3q7CnVTL5sG8fiYwqQZ3YAqwHxiXt9cC95XqDtwAzAJ2NmI9uwKfAnMAAe8Dt7fkfgvQl9eXl+jeyTvHpwIfAN8Ck0vz7mUp4g7c9mHbW/P6MWAX6QLtVjRrAbDG9u+29wN7gdPzLy9OAUn9wJ3Aika4eHdJE0kn+UoA23/Y/okK3DN9wDhJfcB44HsKdXdvReE6uua3nifa/tgps7zBKBSS6+Rue5PtE3nzE6C/NPcu3znAC8DjQHMAsBjvXigigTeRNA2YCWyme9GsKcDBxm4DOdYWL5JOiD8bsRrcLwOOAq/n7p8VkiZQgbvtQ8BzpBfJDgM/295EBe4NenWdkteHxtvmftKdKRTuLmk+cMj2tiFNRXt3o6gELuls4B3gYdu/DPfRDrFWptNImgf8aPuzU92lQ6ytqUB9pEfMV2zPBH5l+Nruxbjn/uIFpMfdi4EJkhYNt0uHWKlTsLq5FncMkpYBJ4BVg6EOHyvCXdJ4YBnwZKfmDrEivIejmAQu6QxS8l5le30OdyuaNUDqxxqkn/T43AbXAfMlHQDWADdJeos63AeAAdub8/Y6UkKvwf1mYL/to7aPA+uBa6nDfZBeXQf4p6uiGW8FSfcC84B7cvcClO1+OekHf1u+XvuBrZIupGzvrhSRwPOo7kpgl+3nG02DRbPg30Wz3gMWShoraTpwBWmgYdSxvdR2v+1pwELgQ9uLqMP9CHBQ0pU5NBf4igrcSV0nsyWNz+fPXNLYSQ3ug/TkmrtZjkmanY95MadYSG6kkXQb8AQw3/ZvjaZi3W3vsH2B7Wn5eh0gTZ44UrL3sLQ9ipp/uK8nPZZsB77Iyx3A+aR/2fZN/nteY59lpJHiPRQyKkwq+jU4C6UKd+AqYEv+7jcA51bk/jSwG9gJvEmaQVCkO7Ca1Fd/nJQ4HvgvrsA1+Xj3AS+TX8ZrwX0vqc948Hp9tTT3Tt5D2g+QZ6GU5N3LEm9iBkEQVEoRXShBEARB70QCD4IgqJRI4EEQBJUSCTwIgqBSIoEHQRBUSiTwIAiCSokEHgRBUCmRwIMgCCrlL+/1s5Nf4bwtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# .. A priori and statistics: \n",
    "#   note: X is temperature, emissivity, y is radiance at nu\n",
    "Xa     = np.array([273., .8])              # a priori X\n",
    "Xfg    = np.array([273., .8])              # first guess X\n",
    "Sa_vec = np.array([3**2., 1.])             # variance for Xa\n",
    "Se_vec = ErrLevelEst**2 * np.ones((Nobs))  # variance in measurement, y\n",
    "\n",
    "# .. Set variables that won't change\n",
    "Sa     = np.diag(Sa_vec)\n",
    "Se     = np.diag(Se_vec)\n",
    "inv_Sa = np.diag(1/Sa_vec);\n",
    "yfg    = greybody(nu,Xfg)\n",
    "yn_1   = yfg + 0.\n",
    "Xn_1   = Xfg + 0.\n",
    "dT     = .1\n",
    "dEps   = .01\n",
    "Kn     = np.ones((Nobs,2))                 # Weighting function\n",
    " \n",
    "# .. Do the retrieval iteratively\n",
    "Niters = 5                          # number of iterations\n",
    "for iter in range(Niters):\n",
    "\n",
    "    # .. Get kernels. Xp is the perturbed Xn_1\n",
    "    #    Temperature, T\n",
    "    Xp      = Xn_1 + 0.; Xp[0] += dT\n",
    "    yp      = greybody(nu,Xp)\n",
    "    Kn[:,0] = (yp - yn_1 ) / dT     # Temperature\n",
    "\n",
    "    #     Emissivity, Eps\n",
    "    Xp = Xn_1 + 0; Xp[1] += dEps\n",
    "    yp = greybody(nu,Xp)\n",
    "    Kn[:,1] = (yp - yn_1 ) / dEps\n",
    "    \n",
    "    # .. Invert as in Rodgers eqn. 5.9 (n-form)\n",
    "    KT_Sem1   = ( solve(Se.T, Kn) ).T\n",
    "    KT_Sem1_K = np.dot( KT_Sem1 , Kn )\n",
    "      \n",
    "    term2 = inv_Sa + KT_Sem1_K \n",
    "    term3 = np.dot( KT_Sem1 , ( yobs-yn_1 + np.dot(Kn,(Xn_1-Xa))  ) )\n",
    "    Xn    = Xa + solve(term2, term3)\n",
    "    \n",
    "    yn    = greybody(nu,Xn)\n",
    "    plt.plot(nu,yn)\n",
    "    #plt.hold(True)\n",
    "    \n",
    "    # .. Set up for next iteration\n",
    "    Xn_1 = Xn + 0.\n",
    "    yn_1 = yn + 0.\n",
    "\n",
    "    print (Xn)\n",
    "\n",
    "plt.plot(nu,yobs,'r-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x123b44d68>,\n",
       " <matplotlib.lines.Line2D at 0x123b44eb8>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZd7/8fd3Jj0kkEAIvYcQegmhiYIUsQHqg6IiqCBgL6uruz5bdH+7j727CiIIUiwggpWmIB0CBIGEJDSRQEKoCSWQcv/+YNxlXZAEMrnnTL6v68p15twzJ+czXMmHk3POnCPGGJRSSjmPy3YApZRSF0cLXCmlHEoLXCmlHEoLXCmlHEoLXCmlHCqgIldWo0YN06hRo4pcpVJKOd66desOGGNifj1eoQXeqFEjkpOTK3KVSinleCLy07nGdReKUko5lBa4Uko5lBa4Uko5lBa4Uko5lBa4Uko5lBa4Uko5lBa4Uko5VIWeB658nykpJv9wLkdys8g/sJeTR7MpOpqDwUWTK24jtm4j2xGVUh5a4JWBMZw6vIejuVnkH9zLycM5FB7NpuTYftwnDhB06iChhYeIKDpCNXOUSCkh8hzfpij9RTaEJVHSbhitew8hODikwt+KUurfpCJv6JCYmGj0k5gV62jOT+ybMooWx9f813MnTDCHpCp57ihOBkZzOqQ6JWE1kPCaBFaNJbRaLSKq16ZazbrkH9pP1uL3aZI1hxgOc4CqpMdeR60rRtG0ZUcL70ypykNE1hljEv9rXAvcTxnDpnkTabjqTwSYIlbWGUFQnVaEVKvtKeU6VI+KJtBdtsMgxUWFpC79jOJ1H9IqfwWBUkxqQAJH42+mZd87qRoV7aU3pFTlpQVeieQdymHHB2Non/c9qe4WBNw0juYt25f7eo7s30PmggnEbv+UBiV7OG6C+bFaH8KSRtCma39cZfzPQSl1blrglcSP331MnR9+T6TJZ3XDMXQZ9gxBQUHeXakxbN/wPUeWT6TFwQWEU8AuqcvPDW+iad+R1KnXyLvrV8rPaYH7ubyjh0j94CG6Hv6Cna6GFA16l7h23Ss8R8Hxo6QtmkrY5unEn95MkXGREtqFkk5307nvEESkwjMp5XRa4H4sZdnXxCx8hNpmP2vr3kH74c8THBJmOxbZO35kz3cTaLRnLjU4zMrIAbQeNY6IyGq2oynlKFrgfij/WD7rP3iCnrkfke2K5cQ1b9Gscz/bsf5LSeFp1k/9Ix13TeAnVz2KbpxIXJsk27GUcozzFbgeZXKo9auXsP/l7lxxYAYbYwcT/fganyxvAFdgEIl3vUTGVR8SafKpN/NaVsx8DVNSYjuaUo6mBe4wx04WMO/dJ2j99Q1UJZ/t/SbR4b4PCAmvajvaBbXofj2ue5exI7QV3Tf/hbWvDiH/6CHbsZRyLC1wB1m3fi27XujJVdnj2Va9N1UeTaZpjxttxyqTqNj6JDy+kNWN7qNT3iIOv9aDbT+utB1LKUfSAneAE6cK+eK9Z2g55xoamL3suOINWj40i5DIGrajXRRXQABd7vw/Mq+eTogpoP6s61n1yYu6S0WpMtIC93ElJYblb43m+qxXyKrakaCHVtOk9wjbscpFi67XEHj/ctJD29M19f+x/pUbyD9y0HYspRxDC9zHLfz0n/TL/4zUBrfR7NFvCYmuZztSuYqKqUPrJ+axqslDtMv/gaOvd2dbylLbsZRyBC1wH5aybhU9Up9he2gbEoa/Dn76IRiX203X4X8j89pPCDBFNJg9mDUf/5/uUlHqArTAfVTuwYNEfjGSU64Qao2cgQR4+ePwPiAhqR/BDyxnS1giSWnPsfGVgeQfybUdSymfpQXug4qKismYcBcNTRbHrhtPeI36tiNVmKgatWj3+DesaPoYrfJXcOz1bmzfsNh2LKV8kha4D1o89e/0OLmE1ISHadBpgO04Fc7ldtH9jr+w7bqZGCM0+PxG1s34G1Tgp4aVcoILFriIxItIyllfeSLyiIhEi8gCEcn0TKMqIrC/W7dsHpfvfI3UiB60ufkvtuNYldD5SkIfXMGPYV3plP4Sayc+qiWu1FkuWODGmHRjTHtjTHugE3ACmA08BSwyxsQBizzz6hLszdpN3YX3ctAdQ5N7PgSX/oEUVT2G9o9/wcpq19H550msmfSE7UhK+YyyNkQfYLsx5idgEDDZMz4ZGFyewSqb06cL2f/BHUSZPMyQyYREVrcdyWe43W6SHpzC6qpXk7T7PVZ9oNsKSkHZC3woMMPzONYYsw/AM61ZnsEqm5UTn6B9YQoZiX+hTkJX23F8jtvtJvHBqayJvIquu95h5eSnbUdSyrpSF7iIBAEDgU/LsgIRGS0iySKSnJurp4Sdy5r5H3FF9iQ2VL+ONtc/aDuOz3IHBNDpoekkR/al2863WPlh5T5GoFRZtsCvBtYbY3I88zkiUhvAM91/roWMMeONMYnGmMSYmJhLS+uHdm9PI375Y+wIaEKrUeNtx/F57oAA2j84g/URvem2/TVWTHvWdiSlrClLgd/Kv3efAMwFfrkoxwhgTnmFqixOnjhOwfRhuMQQPmw6QaHhtiM5QkBgEG0f/JgNVS6ne+bLrJj+D9uRlLKiVAUuImFAP+Czs4afA/qJSKbnuefKP55/S3lvLM2Lt/HT5S8T2yjBdhxHCQgKps1DM0kJv4zuGc+z/KPnbUdSqsKVqsCNMSeMMdWNMUfPGjtojOljjInzTPXK/GWwevbbdDs8lzV1h9P6yttsx3GkgKBgWj08ix/DutFj6z9Y9vFLtiMpVaH0RGMLtm9eTduUv5Ia3JZOd71iO46jBQaFkPDwbDaHJXFZ2t9Y9smrtiMpVWG0wCtY/pGDBM26k2MSTuzd03AHBNqO5HiBwaHEPzyH1NBOdN/yDMtmvmk7klIVQgu8ApmSEra9N4LaJdkcuPpdqsc2sB3JbwQGh9HsoS9ID21P901/Yumsf9qOpJTXaYFXoDUz/h8dji9lbbOHSehS+S5S5W1BoeE0ffhLMkLb0v3HP7J09jjbkZTyKi3wCpK+Zh6dMl5lXVhPut7+Z9tx/FZQaBUaP/QF20Na0y3lKZbOed92JKW8Rgu8AhzZv4for8ewzxVLs3smI3qRKq8KDqtKw4e+ZGdIC7quf4KlcyfZjqSUV2iTVIAd0x4hwhzj5A2TqRqlF6mqCMHh1aj/4NfsDo6jy7rfsfTLKbYjKVXutMC9bE/6etofWci6WkNp3raL7TiVSkiVKOo++A17gpuStPZRVn//he1ISpUrLXAvO/DlXzlBCPE36dXzbAiJiKbOA99wICCW+MVjSduSYjuSUuVGC9yLdqeuon3+ElLq3UqNmrVtx6m0QiJrEDJiFiIQ9umtZO3Lsh1JqXKhBe5FR756hqMmnFY3/tF2lEqveoME8gZNpjb7OTDhFvKOH7cdSalLpgXuJTs3/kDb4yvY1GAYUdX1Mrq+oH6Hvuzs/jztijex4Z93UVhUbDuSUpdEC9xLjn/7LIdNBG1u1Nt/+ZL4/qPYHDeWK47PY9GEP2D0JsnKwbTAvWDHuoW0PrmWzY3vpGpUtO046lda3/YcqdX7MyB7HItm6U00lHNpgXvBqQV/4yBVaXfj47ajqHMRocWYKWwPacVlm55m5ZJ5thMpdVG0wMtZ5uqvSChIIbXpKCIjq9mOo87DFRRK3bGzOeKuTtx395Causl2JKXKTAu8PBlDyaK/k0M0HW98zHYadQEh1WIJHjGTYCkk+JNbycrOth1JqTLRAi9H6Ss+J/70FjKajyE8vIrtOKoUohq24ej1E2nAXvZNuJW8EydtR1Kq1LTAy4sxuBf/g73EkHjDw7bTqDKo1+lqfuryLIlF61nz9ig9vVA5hhZ4OUlb8gnNCjPYnnAfoaGhtuOoMmp29QOkNbmLvse/ZN7Ev+rphcoRtMDLgSkpJmTZc/xMLToPus92HHWREoa9Qnp0L67JepN5s/UStMr3aYGXg7TvptO4aAc/tXmQkJAQ23HUxXK5iBsznZ9DmtNz41MsX7rIdiKlfpMW+CUyxUVUWfkCu6QuSdePsR1HXSJXcDixY2Zz0h1B04Wj2LJ1q+1ISp2XFvgl2rLgAxoU7yar3SMEBekd5v1BSHRd3Hd8SqScwP3RULJycm1HUuqcSlXgIlJNRGaKyFYRSRORbiISLSILRCTTM43ydlhfY4oLiVrzMtulIUnX3W07jipHUY07cuTaccSxi5/eu428EwW2Iyn1X0q7Bf468K0xpgXQDkgDngIWGWPigEWe+Upl8zfvUbdkL/s7/Y7AgADbcVQ5q9N5MLs6/5nuRWtY9s/7KCousR1Jqf9wwQIXkUjgcuB9AGPMaWPMEWAQMNnzssnAYG+F9EUlhaepse410l1N6TxgmO04ykuaXvsYGY1u45pjs/hmygu24yj1H0qzBd4EyAUmicgGEZkgIuFArDFmH4BnWvNcC4vIaBFJFpHk3Fz/2Ze4+et/UtvkcCjpcQIC3LbjKC9qfsebbIvswlW7XmDJgjm24yj1L6Up8ACgI/COMaYDcJwy7C4xxow3xiQaYxJjYvzjxgbFpwuolfImqe54kvoNtR1HeZs7gIZjPuZAQC3aLLuftDS98JXyDaUp8D3AHmPMas/8TM4Ueo6I1AbwTPd7J6Lv2fzF69Q0B8jv9iRut57IUxkEhkcROmImAVJC4Ce3k3vwoO1ISl24wI0x2cDPIhLvGeoDpAJzgRGesRFApfjbsqjgGPU2v8OmgNZ07n2D7TiqAkU1aMmhq8fRuGQ3u8bfTmFRke1IqpIr7ebjg8A0EfkRaA/8A3gO6CcimUA/z7zf2zznVaqbwxT0/AMu3fqudBp1uZ7Utk/R+dRKVrynlwxWdpXq3DdjTAqQeI6n+pRvHN92+kQeDdPGkxLYnsTLr7UdR1nS5sYn2ZC9hStyJrPy85Z0GzzWdiRVSekmZBmkfv4iUeRR3OuPiIjtOMoWEdrc8x5bg1vTYcP/kpa82HYiVUlpgZfSqWOHaZIxkeSgznTs3t92HGVZQFAIte/5lMOualT/8i5y9+6yHUlVQlrgpZQ2+3kiOYa7z//q1rcCoGqNOhT8zzSqmOMcmTiEUyeP2Y6kKhkt8FIoyDtAs+2TWR3cnfZJV9iOo3xI41ZdSOv2MnFFGaS+eyemRD9uryqOFngppM1+gTBzkpD+uvWt/lunAXewtP5YOhxdwIaP/mo7jqpEtMAvwBQXUn/Xp6QEd6Jdpx624ygf1f3O/2N1eG/ap79BxpKPbcdRlYQW+AVkLPuMGuYQBW2H246ifJjb7aLFmClkuJtS9/uH2L9tve1IqhLQAr+AwjWTyKUaHfrcYjuK8nFVIyMJHjaD4yYUM30oBUcrzdUllCVa4L/hSPYuEo6tIjV2IKGheq9LdWGNmzRnZ9/xVCs+xJ5xQzBFp2xHUn5MC/w3bJ//Lm4x1Ok92nYU5SBdevbn+/i/0OxECumT7gVjbEdSfkoL/DxMcRH1ds4kJbA9cS3a2I6jHKb/0Af4ptqttMiaxY6vX7UdR/kpLfDzyFw5l1iTy/E2ercdVXYul9Bz7BusDEiiwdq/sT/lW9uRlB/SAj+PgtUTOWgiad/3dttRlENVCQmi9t1T2UldQufczcmcTNuRlJ/RAj+HvNyfSchbwZaa1xIeFmY7jnKwRnViOXD9ZIpL4ODEWzCnj9uOpPyIFvg5bJs/nkApplavMbajKD/QrVMnlrb5P+oU7GDbxHv0oKYqN1rgv2JKiqm1/RM2BbSmeasOtuMoP3HdTcP5qvoI4rK/YtvXr9uOo/yEFvivbF/7DXVKsslrqfu+VfkREXqPfpHVAYk0XPssOVt+sB1J+QEt8F85vvx9jphw2va/w3YU5WeqhARR684pZFMD98wRnDy0z3Yk5XBa4Gc5dmgfCUeXsLnGNURUibAdR/mhhvXqsu+qcVQpyWfPhKGY4kLbkZSDaYGfJX3+ewRJMdUvv8d2FOXHkrr1Zmn808SdSGHTlN/ZjqMcTAv8F8ZQM+NjtrgTaNE2yXYa5ef6DH2ExZEDafvTZLZ+N9V2HOVQWuAe25MXUL9kD4db3Ko3bVBe53IJiWPHkeqOp/4PvyN7e4rtSMqBtMA98lZMIN+E0qb/CNtRVCVRJSyMsGFTKTBBFE6/nYJjR2xHUg5TqgIXkV0isklEUkQk2TMWLSILRCTTM43yblTvOX4kl5aHvmNj9FVUrVrNdhxViTRq3Jyfer9FnaIsto4brvfUVGVSli3w3saY9saYRM/8U8AiY0wcsMgz70jp8ycQLIVEXaYHL1XF69hrECsaP0j7/CWsnfE323GUg1zKLpRBwGTP48nA4EuPY4ExVE+fwVZXHC076j0vlR097niGdeE96ZjxGqkrvrIdRzlEaQvcAPNFZJ2I/HJ3g1hjzD4Az7TmuRYUkdEikiwiybm5uZeeuJztSllMw+KfOBA/VA9eKmtcbhfNR08hy12H2Pn3kr1nh+1IygFKW+A9jDEdgauB+0Xk8tKuwBgz3hiTaIxJjImJuaiQ3nRo6XscN8G06n+37SiqkouoGg03TyXEnOLwB7dSUHDSdiTl40pV4MaYvZ7pfmA2kATkiEhtAM/UcXdwPZl3mBaHFrKxWj+ioqJtx1GKhi06sK378yQUbSV53L0YvXKh+g0XLHARCReRiF8eA/2BzcBc4Jdz7kYAc7wV0lu2LnifME4R0WOk7ShK/Uu7q+4kuc7tXHZ4Nss/e8t2HOXDSrMFHgssE5GNwBrgK2PMt8BzQD8RyQT6eeYdpWradLa5GtM6sZftKEr9h453v056SFs6/fgsm5KX2Y6jfNQFC9wYs8MY087z1coY83fP+EFjTB9jTJxnesj7ccvPT5uX06RoO9nNbkFc+nkm5VtcAYHUuecjjruqEPXl3WTn6JUL1X+rtM11YPE4TpogWl41ynYUpc4ponpdCm6YRKw5wO4Jwyk4rVcuVP+pUhZ4wfGjtDgwj41VexNd3ffOjFHqF/Xa9iKz49MkFa5h6bhH9KCm+g+VssDTFkwmnAJCu+rBS+X7Wg58jM21BtPv4FS++/Rt23GUD6mUBV5ly1R2Sn3adOlnO4pSFyZCq1HjyQxty2Vb/kry8vm2EykfUekKfE/aGuIK08lqcjMud6V7+8qhJCCYemNmcdgdTcP597Bze7rtSMoHVLoGy/5+HKdMIPF68FI5TGi1mrhu/5gwKeD0tKEcOXLYdiRlWaUq8FMn84nf/zUbIy4npmYd23GUKrOaTTuwr+/bxBXvJH3cHRQWFdmOpCyqVAWeuvBDIjhBUNJdtqModdGaXfY/bG71O7qcXMqyCU/YjqMsqlQFHrZpGrulNm16XGs7ilKXpO2Q/2Vj9WvonT2RZZ+Ptx1HWVJpCnxv5gbiT29md6MhuPXgpXI6EVqPmURmUEs6bXiaTWsX206kLKg0TZa1aBynjZvm/Udf+MVKOYA7KITY0bPIc1Wl5ld3kfXzTtuRVAWrFAVeeOoEcdlfsrFKD2rWrm87jlLlJrJGHYpumU4Exzn2wRCOHcu3HUlVoEpR4KmLplGNfAIS9eCl8j91WySx6/LXiC/OZMs7wykp1hsjVxaVosBdmz5mLzG06TnQdhSlvKLllbeR3PRBuhz/juUf/NF2HFVB/L7A8w5kk3BiHTtrDSAgIMB2HKW8ptPtz7K+aj96/vwOa76ZfOEFlOP5fYFnLJ5KgJQQ0+0221GU8ipxuWh97xQyA+NpveoJMjYutx1JeZnfF3hYxhx2ST3i2nS1HUUprwsKCSN65EzypQqRs4eTu2+37UjKi/y6wA9k7aTFqU3srXeN3nVHVRrVazXg+I0fUtXkcfD9myk4ecJ2JOUlft1q2xd/iEsMdXsOsx1FqQrVpG0PtnZ7gRZFaWx6905MiZ6Z4o/8usCjdn5BprspDZu3sx1FqQrXYcBdrGowms5H57F62l9tx1Fe4LcFvnf7FpoXZXCg0fW2oyhlTZc7n2ddlV4kbXuDDfOm2I6jypnfFvjupVMBaNzrDstJlLJHXC5a3juVjMB4ElY8RsbahbYjqXLktwVea/eXpAW2olb9ZrajKGVVaHgEMaM/44CrBjW/upOsbT/ajqTKSakLXETcIrJBRL70zEeLyAIRyfRMo7wXs2x2pK6hUclu8poNsh1FKZ9QvWZdSm77lBKAaUM4vD/LdiRVDsqyBf4wkHbW/FPAImNMHLDIM+8TcpZNo8i4aNbrdttRlPIZDeLakHPtZKJLDnFg/GAKjufZjqQuUakKXETqAdcCE84aHgT88nndycDg8o12cUxJCQ32fUNaaAeqx9azHUcpn5LQuQ9bur9Kk8JMMt6+mZKiQtuR1CUo7Rb4a8DvgbNPJo01xuwD8ExrnmtBERktIskikpybm3tJYUsjfcMS6pocCuJv8Pq6lHKixKuGsSr+SdqeWMn6caPBGNuR1EW6YIGLyHXAfmPMuotZgTFmvDEm0RiTGBMTczHfokwOr57BaRNAi963en1dSjlV91ufYnnN20nM/YzkaX+xHUddpNJsgfcABorILuAj4EoRmQrkiEhtAM90v9dSllJRYSHN9s8ntUoXIqrVsB1HKZ8lInQd8yZrwnuTuO11Nn79nu1I6iJcsMCNMX8wxtQzxjQChgLfGWOGAXOBEZ6XjQDmeC1lKaWtnkcMhzGtbrIdRSmf53a7aXP/dDYHtiFh9ZNkrP7adiRVRpdyHvhzQD8RyQT6eeatOr7uY06YYBKuGGI7ilKOEBoWRp0xs8hy1aH2NyPJSl9vO5IqgzIVuDFmsTHmOs/jg8aYPsaYOM/0kHcils6pUwXEH/6erVUvIyQ80mYUpRwlukYsAXfMooAgAj4awuEcvQStU/jNJzG3LJ1DFPkEttetb6XKqn6TeHIHTqVKyTGOvDeIgmOHbUdSpeA3BV648VPyCKdFD584HV0px2nZsSdpPd+kfuEudrx9E8WFp21HUhfgFwV+/Fg+rfOWkhHdm8DgUNtxlHKsxL43s7Lln2l5ch2b3r1TzxH3cX5R4FuWzCRcCqiSONR2FKUcr+ctj7Kk9kjaH/yK9VOetB1H/Qa/KHDXllkcpBrNk662HUUpv3DZqJdYHjGAjjvHsemLt2zHUefh+AI/fOggrY+vYmdsP1wBAbbjKOUX3G4Xne6fzIagjiQk/4nMFbNtR1Ln4PgCT1v8ESFSSFQX/ei8UuUpJCSEhmNnssvdgDrzx/Lzph9sR1K/4vgCD02fTbbUpEn73rajKOV3oqOrEzz8Mw5TlchZt7Jv61rbkdRZHF3gOdlZtC5Yz546AxCXo9+KUj6rfqOmnLrtc04STMhHN5K7Y6PtSMrD0a2X+f00AqWYWj30xg1KeVPT5i05OmQWxUaQDwdz8OettiMpHF7gkTvmssdVj3oJXWxHUcrvxbfqQPbgj3GXFFI08XqOZu+wHanSc2yB7961jdanN5PT8DoQsR1HqUqhdYdu7LpmKqElxzg2/lqOHdhjO1Kl5tgC3/XDNFxiqH/5MNtRlKpUOnTpRXqfSVQrPsihd6/h5BHrtwKotBxZ4MYYYnZ9yY6AptRs3MZ2HKUqnc6XD2Bjz3HULNxLzttXc+qY1YuRVlqOLPBt6ZtIKMngSNOBtqMoVWl173sDa7q8QZ3TO9nz5jUUnjhqO1Kl48gCz1o2DYAmvYZbTqJU5Xb5NbexvMOLNCxIZ9ebAyk+dcJ2pErFcQVeUmKon/U1GcGtqFa7ie04SlV6vQePZHGrv9H0xEa2vTmYktMFtiNVGo4r8NSNq2hqdnOy+SDbUZRSHn1vfoCFzf5I/LHVpL99M6ZIryVeERxX4AdXzaDYCM1632E7ilLqLP2GPcG39R8j4egS0t65HUqKbUfye44q8KKiYprkzCMjvBPh0XVsx1FKnUVEuOruP/NNrbG0PDif1PF3QUmJ7Vh+zVEF/uPaxdQnm6KEG2xHUUqdg4jQf/RzfFt9BC2z55A66T69q48XOarA85M/ohA3cb1usx1FKXUebpfQ995XmV91CC1/nkHa1Me0xL3EMQVecLqQFgcWkBHRlZCIaNtxlFK/ISDAzRUPvMuiKteTsH0iaR//yXYkv3TBAheREBFZIyIbRWSLiDzjGY8WkQUikumZRnkz6MblXxMrh3G3HeLN1SilyklwYADdH5zE4tB+JGx9k60f/VG3xMtZabbATwFXGmPaAe2BASLSFXgKWGSMiQMWeea95lTKp5wkmGaX/Y83V6OUKkehwYEkPjSNxaH9aLH1bVKnPKolXo4uWODmjGOe2UDPlwEGAZM945OBwV5JCOQfP0GbI4vZVq0nAaER3lqNUsoLqoQG0/XRGSyKGEjLnZPYMmG0np1STkq1D1xE3CKSAuwHFhhjVgOxxph9AJ5pzfMsO1pEkkUkOTc396JC/rh0DlGST2inWy5qeaWUXSFBgfR86APmRw2lVdYnbHlnGKa40HYsxytVgRtjio0x7YF6QJKItC7tCowx440xicaYxJiYmIsKGZo+h3zCadpVL16llFMFBbrp88A7zKs5ila5X5H61s2UFJ6yHcvRynQWijHmCLAYGADkiEhtAM/UaxcF7jjqbbhlGhIY4q1VKKUqgNvtov+9LzGv3kO0Ovwd6W8Movj0SduxHKs0Z6HEiEg1z+NQoC+wFZgLjPC8bAQwx1shCa9ORILedV4pfyAi9B/5LAub/oH4vFVse/VqTp/Isx3LkUqzBV4b+F5EfgTWcmYf+JfAc0A/EckE+nnmlVLqgkSEvnc8xQ+eqxj+9PpVnMzTm0KUlZgKPKUnMTHRJCcnV9j6lFK+b+nciXRZ9zh7AhtT876vqBJdy3YknyMi64wxib8ed8wnMZVS/qnnwLtZ3/2f1Cn8iUNv9+NIzm7bkRxDC1wpZV3Xq4aSeuVEqhflcPzdfuTuybQdyRG0wJVSPqHjFQPZfvU0qpTkU/L+APbt2Gw7ks/TAldK+Yy2XfuRfcOnBJpTBE65lt1pa21H8mla4EopnxLfvgdHbpmDQYj4+Aa2pSy1HclnaYErpXxOk4ROFNzxFQWEEDt7CGmr59uO5JO0wJVSPql+01bIyG854hn4FqAAAAoeSURBVIqi4dfDWL9guu1IPkcLXCnls2rVb0bY2PlkBTag/bL7WDHlzxi9kuG/aIErpXxa9dj61H/se1Iie9F9x+usff1WCk6esB3LJ2iBK6V8XkhYBB0e/YzVDe4h6ei37HilLwdysmzHsk4LXCnlCOJy0eXul9iQ9ApNTmdw+t1e7NiyxnYsq7TAlVKO0uGakewZ/BmBppDYT65nw8KPbEeyRgtcKeU4zTpcDvd8z76AurRbOpYVHz5TKQ9uaoErpRwppm5j6j66mJSIy+m+/RXWvHE7p05VrptDaIErpRwrtEokHR6bzar6I+ly5Gu2vdSXg/v32o5VYbTAlVKOJi43XUe+wvrOL9LsdDoF7/RiR2rluO+AFrhSyi90vHY0uwd+Sog5RczH17Fh0Se2I3mdFrhSym/EdepN8ajv2B9Qm7Y/jGbFtGf9+uCmFrhSyq/UrNeUOo8uYVOVHnTPfJnVbw7324ObWuBKKb8TWiWSto/NZXW9u+h6+AsyX+5P7j7/u1WbFrhSyi+53G66jHqN9Z2eJ+5UGjKuJxsXz7Idq1xpgSul/FrH68eSfcs3HHNF0m7x3ax6ZyynCvzjYlha4Eopv9ewZWdqPb6S1TVupGvODH5+8TL2ZG60HeuSXbDARaS+iHwvImkiskVEHvaMR4vIAhHJ9EyjvB9XKaUuTkhYFbo8MIn13d+mRnEO0VP7se7zN8AY29EuWmm2wIuA3xljEoCuwP0i0hJ4ClhkjIkDFnnmlVLKp3XsP4yCUUvZEdyCTil/YsMrN3Ds6EHbsS7KBQvcGLPPGLPe8zgfSAPqAoOAyZ6XTQYGeyukUkqVp1r1mpDw++9Y3vB+2uQtIf+1Lmxbt9B2rDIr0z5wEWkEdABWA7HGmH1wpuSBmudZZrSIJItIcm5u7qWlVUqpcuIOCKDHXf8g/dqZlBgXjeYOIXnyk5QUFdmOVmqlLnARqQLMAh4xxuSVdjljzHhjTKIxJjEmJuZiMiqllNe0SupD+MMrWRfRm8Sd75L5Yi8OZm23HatUSlXgIhLImfKeZoz5zDOcIyK1Pc/XBvZ7J6JSSnlXtajqJP1uFsvb/J16BZkEvteT1IUf2o51QaU5C0WA94E0Y8wrZz01FxjheTwCmFP+8ZRSqmKICD1ueoCc2xaw112HlsseYMPbwzl9It92tPMqzRZ4D+AO4EoRSfF8XQM8B/QTkUygn2deKaUcrUl8Wxo+vpTFNYfRbv9ccl7uxt6tvnnvTTEVeA5kYmKiSU6uHNfpVUo535pFn9Fo6WNUM/mkxo2m7c1/wRUUUuE5RGSdMSbx1+P6SUyllDqPpD43UjxmGclhl9F+2z/Z93xHfl73re1Y/6IFrpRSv6F27Xp0+/3nLOs6jpLiYup/cQtb3rqFgsP7bEfTAldKqQsRES4bMJTwR9awIGYEcbkLKHy9E5lfvwEWbxihBa6UUqUUXa0q/e5/gy2DvmWbuzFxa/7Ezhd6cHjHOit5tMCVUqqMOnRMIuHJJXwT9wyRJ/cQMaUvWyc/hDlVsaccaoErpdRFCAkK4OrbH+HoyOV8H3oVLXZO5uDz7dm3amaFZdACV0qpS9CkQQP6PDGDhd0+5FBJGLW/Hcm2N66n4MAur69bC1wppS6RyyX0vWog0Y+uZE7MWOocXI15K4ldc/8BxYXeW6/XvrNSSlUyNapWYdD9z7Np8ALWudvRaP3z7H0hibyMZV5Znxa4UkqVsy4d2tHpyW+Z1fwFKDhK5PRrSf9+ermvRwtcKaW8IDTIzU23jSF/1HJmRg6napsB5b6OgHL/jkoppf4lvn5t4h970yvfW7fAlVLKobTAlVLKobTAlVLKobTAlVLKobTAlVLKobTAlVLKobTAlVLKobTAlVLKoSr0psYikgv8dJGL1wAOlGOciqTZ7XBqdqfmBs3uLQ2NMTG/HqzQAr8UIpJ8rrsyO4Fmt8Op2Z2aGzR7RdNdKEop5VBa4Eop5VBOKvDxtgNcAs1uh1OzOzU3aPYK5Zh94Eoppf6Tk7bAlVJKnUULXCmlHMonClxE6ovI9yKSJiJbRORhz3i0iCwQkUzPNOqsZf4gIttEJF1ErrKX/l953CKyQUS+9Mw7IruIVBORmSKy1fPv381B2R/1/LxsFpEZIhLiq9lFZKKI7BeRzWeNlTmriHQSkU2e594QEbGU/UXPz8yPIjJbRKr5WvZz5T7rucdFxIhIDV/LXSbGGOtfQG2go+dxBJABtAReAJ7yjD8FPO953BLYCAQDjYHtgNvye3gMmA586Zl3RHZgMjDK8zgIqOaE7EBdYCcQ6pn/BLjTV7MDlwMdgc1njZU5K7AG6AYI8A1wtaXs/YEAz+PnfTH7uXJ7xusD8zjzocIavpa7LF8+sQVujNlnjFnveZwPpHHmF3QQZwoGz3Sw5/Eg4CNjzCljzE5gG5BUsan/TUTqAdcCE84a9vnsIhLJmR/y9wGMMaeNMUdwQHaPACBURAKAMGAvPprdGPMDcOhXw2XKKiK1gUhjzEpzplmmnLVMhWY3xsw3xhR5ZlcB9Xwt+3n+zQFeBX4PnH0Gh8/kLgufKPCziUgjoAOwGog1xuyDMyUP1PS8rC7w81mL7fGM2fIaZ34gSs4ac0L2JkAuMMmz+2eCiITjgOzGmCzgJWA3sA84aoyZjwOyn6WsWet6Hv963La7ObNlCj6eXUQGAlnGmI2/esqnc5+PTxW4iFQBZgGPGGPyfuul5xizcj6kiFwH7DfGrCvtIucYs3UuZwBn/sR8xxjTATjOmT/lz8dnsnv2Fw/izJ+7dYBwERn2W4ucY8xXz6E9X1afew8i8jRQBEz7ZegcL/OJ7CISBjwN/PlcT59jzCdy/xafKXARCeRMeU8zxnzmGc7x/AmDZ7rfM76HM/uxflGPM38+29ADGCgiu4CPgCtFZCrOyL4H2GOMWe2Zn8mZQndC9r7ATmNMrjGmEPgM6I4zsv+irFn38O9dFWePWyEiI4DrgNs9uxfAt7M35cx/+Bs9v6/1gPUiUgvfzn1ePlHgnqO67wNpxphXznpqLjDC83gEMOes8aEiEiwijYE4zhxoqHDGmD8YY+oZYxoBQ4HvjDHDcEb2bOBnEYn3DPUBUnFAds7sOukqImGen58+nDl24oTsvyhTVs9ulnwR6ep5z8PPWqZCicgA4ElgoDHmxFlP+Wx2Y8wmY0xNY0wjz+/rHs6cPJHty7l/k+2jqJ7/uC/jzJ8lPwIpnq9rgOrAIiDTM40+a5mnOXOkOB0fOSoM9OLfZ6E4IjvQHkj2/Nt/DkQ5KPszwFZgM/AhZ84g8MnswAzO7Ksv5ExxjLyYrECi5/1uB97C82lqC9m3cWaf8S+/r+/6WvZz5f7V87vwnIXiS7nL8qUfpVdKKYfyiV0oSimlyk4LXCmlHEoLXCmlHEoLXCmlHEoLXCmlHEoLXCmlHEoLXCmlHOr/A20TydS3JoqEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(nu,yn,nu,yobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD5CAYAAAAqaDI/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAbyUlEQVR4nO3df5Dcd33f8ecrJwnONu7ZGM7SSROJVqhRS2tZF2Fwm54xRD/CWBpKO3IKCDdUQ8ZigCQiUjxNJn90EIhJDBMPGo1xYgca1QFFqHDhMLa2bUrt2kZgWRYXHYJadyewcTjDocP64Xf/2O/Zq/Pe3e5+98d39/t6zOzcfj/fz3f3Jc3tvm8/389+vooIzMwsv36p1QHMzKy1XAjMzHLOhcDMLOdcCMzMcs6FwMws51wIzMxybkE9HkTSBuDTQBdwd0TsmbFfyf5NwFng/RHxrWTfR4EPAAEcA26LiF/M9XzXXHNNLF++vKasP//5z7n88strOrbVnL012jV7u+YGZ2+Uxx9//McR8bpX7IiIVDeKb/7fA94ALAK+A6ye0WcT8LeAgBuAR5L2PuD7QHeyfT/FIjHnc65duzZqdeTIkZqPbTVnb412zd6uuSOcvVGAx6LMe2o9hobWASMRcSoizgEHgM0z+mwG7kuyPAz0SFqc7FsAdEtaAFwGjNchk5mZVagehaAPOF2yPZq0zdsnIsaATwFPA2eA5yPi63XIZGZmFarHOQKVaZu5bkXZPpKuovhpYQUwAfy1pPdExOdf8STSdmA7QG9vL4VCoaawk5OTNR/bas7eGu2avV1zg7M3Wz0KwSiwrGR7Ka8c3pmtz9uB70fEswCSDgJvBV5RCCJiP7AfoL+/PwYGBmoKWygUqPXYVnP21mjX7O2aG5y92eoxNPQosFLSCkmLgK3A4Rl9DgPvU9ENFIeAzlAcErpB0mXJzKKbgRN1yGRmZhVK/YkgIi5I2gEMUZxBdE9EHJf0wWT/PmCQ4syhEYrTR29L9j0i6YvAt4ALwFGSv/rr7dDRMfYODTM2MUXfww+xc/0qtqyZeSrDzCx/6vI9gogYpPhmX9q2r+R+ALfPcuwfAX9UjxyzOXR0jN0HjzF1/iIAYxNT7D54DMDFwMxyLxffLN47NPxSEZg2df4ie4eGW5TIzCw7clEIxiemqmo3M8uTXBSCJT3dVbWbmeVJLgrBzvWr6F7YdUlb98Iudq5f1aJEZmbZUZeTxVk3fUL4pVlDPd2eNWRmlshFIYBiMdiypq8tv+xh7W96+vL4xBRL/IeIZUxuCoFZq3j6smVdLs4RmLWSpy9b1rkQmDWYpy9b1rkQmDWYpy9b1rkQmDWYpy9b1vlksVmDlU5f9qwhyyIXArMmmJ6+bJZFHhoyM8s5FwIzs5xzITAzy7m6FAJJGyQNSxqRtKvMfkn6TLL/CUnXl+zrkfRFSd+VdELSW+qRyczMKpO6EEjqAu4CNgKrgVslrZ7RbSOwMrltBz5bsu/TwNci4p8C/xJfs9jMrKnq8YlgHTASEaci4hxwANg8o89m4L4oehjokbRY0pXArwGfA4iIcxExUYdMZmZWoXoUgj7gdMn2aNJWSZ83AM8Cfy7pqKS7JV1eh0xmZlahenyPQGXaosI+C4DrgQ9FxCOSPg3sAv7zK55E2k5xWIne3l4KhUJNYScnJ2s+ttWcvTXaNXu75gZnb7Z6FIJRYFnJ9lJgvMI+AYxGxCNJ+xcpFoJXiIj9wH6A/v7+qPWaAu18PQJnb412zd6uucHZm60eheBRYKWkFcAYsBX4zRl9DgM7JB0A3gw8HxFnACSdlrQqIoaBm4Gn6pDJZvCFUcxsNqkLQURckLQDGAK6gHsi4rikDyb79wGDwCZgBDgL3FbyEB8CviBpEXBqxj6rA18YxczmUpe1hiJikOKbfWnbvpL7Adw+y7HfBvrrkcPKm+vCKC4EZuZvFueAL4xiZnNxIcgBXxjFzObiQpADvjCKmc3F1yPIAV8Yxczm4kKQE74wipnNxkNDZmY550JgZpZzLgRmZjnnQmBmlnMuBGZmOedCYGaWc54+apYDXn3W5uJCYNbhvPqszcdDQ2Ydbq7VZ83AhcCs43n1WZuPC4FZh/PqszafuhQCSRskDUsakfSKaw6r6DPJ/ickXT9jf5eko5K+Uo88ZvYyrz5r80ldCCR1AXcBG4HVwK2SVs/othFYmdy2A5+dsf/DwIm0Wczslbas6ePj73oTfT3dCOjr6ebj73qTTxTbS+oxa2gdMBIRpwCSC9Rv5tKL0G8G7ksuWfmwpB5JiyPijKSlwG8A/wX4nTrkMbMZvPqszaUehaAPOF2yPQq8uYI+fcAZ4E7gY8Br5noSSdspfpqgt7eXQqFQU9jJycmaj201Z2+Nds3errnB2ZutHoVAZdqikj6S3gk8ExGPSxqY60kiYj+wH6C/vz8GBubsPqtCoUCtx7aas7dGu2Zv19zg7M1Wj5PFo8Cyku2lwHiFfW4EbpH0A+AA8DZJn69DJjMzq1A9CsGjwEpJKyQtArYCh2f0OQy8L5k9dAPwfESciYjdEbE0IpYnxz0UEe+pQyYzM6tQ6qGhiLggaQcwBHQB90TEcUkfTPbvAwaBTcAIcBa4Le3zmjWT1+qxTlaXtYYiYpDim31p276S+wHcPs9jFIBCPfKY1dM3x8/zlw96rR7rXP5msdk8vvT3571Wj3U0FwKzeTz3i5mT4Iq8Vo91ChcCs3m89tXlZj97rR7rHC4EZvP4t29c6LV6rKO5EJjN461LFnqtHutovkKZWQW8Vo91Mn8iMDPLORcCM7OccyEwM8s5FwIzs5xzITAzyzkXAjOznHMhMDPLORcCM7OccyEwM8u5uhQCSRskDUsakbSrzH5J+kyy/wlJ1yftyyQdkXRC0nFJH65HHjMzq1zqQiCpC7gL2AisBm6VtHpGt43AyuS2Hfhs0n4B+N2I+BXgBuD2MseamVkD1eMTwTpgJCJORcQ5iheh3zyjz2bgvih6GOiRtDi5bvG3ACLiZ8AJwAu6mJk1UT0KQR9wumR7lFe+mc/bR9JyYA3wSB0ymZlZheqx+mi5q3bMvKTTnH0kXQF8CfhIRPy07JNI2ykOK9Hb20uhUKgp7OTkZM3Htpqzt0a7Zm/X3ODszVaPQjAKLCvZXgqMV9pH0kKKReALEXFwtieJiP3AfoD+/v4YGBioKWyhUKDWY1vN2VujXbO3a25w9marx9DQo8BKSSskLQK2Aodn9DkMvC+ZPXQD8HxEnJEk4HPAiYj4kzpkMTOzKqX+RBARFyTtAIaALuCeiDgu6YPJ/n3AILAJGAHOArclh98IvBc4JunbSdsfRMRg2lxmZlaZulyhLHnjHpzRtq/kfgC3lznu7yh//sDMzJrEl6q0tnDo6Bh7h4YZn5hiSU83O9ev8qUjzerEhcAy79DRMXYfPMbU+YsAjE1MsfvgMQAXA7M68FpDlnl7h4ZfKgLTps5fZO/QcIsSmXUWFwLLvPGJqarazaw6LgSWeUt6uqtqN7PquBBY5u1cv4ruhV2XtHUv7GLn+lUtSmTWWXyy2DJv+oSwZw2ZNYYLgbWFLWv6/MZv1iAeGjIzyzkXAjOznHMhMDPLOZ8jMDObR6cvceJCYE3xzfHz3LHnoY59IVnnysMSJy4E1nCHjo7xF0+e49yLxe1OfCFZtqX5i36uJU465ffX5wis4fYODb9UBKZ5rSBrlum/6Mcmpghe/kPk0NGxio7PwxInLgTWcHl4IVl2pV20MA9LnNSlEEjaIGlY0oikXWX2S9Jnkv1PSLq+0mOt/eXhhWTZlfYPkTwscZK6EEjqAu4CNgKrgVslrZ7RbSOwMrltBz5bxbHW5nauX8WiGb9pnfZCsuxK+4fIljV9fPxdb6KvpxsBfT3dfPxdb+qY8wNQn5PF64CRiDgFIOkAsBl4qqTPZuC+5JKVD0vqkbQYWF7Bsdbmtqzp46kTT/HVp7s8a8iabuf6VZfM+oHq/xDp9CVO6lEI+oDTJdujwJsr6NNX4bF185GPfIRCoUBPT0+jnqKhJiYm2j77Lyfbd34N7mxposq16/97PXP/ePIFTv/DFC9cuMirFnSx7OpurrniVXV57HLq/X/eNfkCPynJf8XV3dz5v17VkN/BRv++XHfdddx5Z32T16MQlLv4fFTYp5Jjiw8gbac4rERvby+FQqGKiEWjo6NcvHiRiYmJqo/NAmdvjXbNXq/cPz0X/PDnL/Ji8sp84cJFTj07ydmzZ7lyUbmXcHr1/j9fAKy4EiAZ678wxUSDJis0+vdldHS0pve/udSjEIwCy0q2lwLjFfZZVMGxAETEfmA/QH9/fwwMDFQddGBggEKhQC3HZoGzt0a7Zq9X7hv3PMSLZd40F/d08793vS3145fTrv/n0J7Z6zFr6FFgpaQVkhYBW4HDM/ocBt6XzB66AXg+Is5UeKyZtZCn/3a+1J8IIuKCpB3AEMXPXfdExHFJH0z27wMGgU3ACHAWuG2uY9NmMrP6WdLTzViZN31P/+0cdVliIiIGKb7Zl7btK7kfwO2VHmtm2VGPWTeWbV5ryMzm5EuFdj4XAjObV6fPo887rzVkZpZzLgRmZjnnQmBmlnMuBGZmOeeTxZYLnX7NWbM0XAis4+XhmrNmaXhoyDpe2itUmXU6FwLreF4rx2xuLgRt4tDRMX63cJYVu77KjXseqvjC2+ZLZZrNx4WgDUyPcT/3iyB4eYzbxaAyebjmrFkaLgRtwGPc6eThmrNmaXjWUBvwGHd6XivHbHb+RNAGPMZtZo3kQtAGPMZtZo2UqhBIulrSA5JOJj+vmqXfBknDkkYk7Spp3yvpu5KekPQ3knrS5OlU02Pcr321PMZtZnWX9hzBLuDBiNiTvMHvAn6/tIOkLuAu4B0UL2L/qKTDEfEU8ACwO7lk5SeA3TOPt6Ita/roef5k210U28yyL20h2AwMJPfvBQq88o18HTASEacAJB1IjnsqIr5e0u9h4N0p8zSM16oxs06l4uWEazxYmoiInpLtn0TEVTP6vBvYEBEfSLbfC7w5InbM6Pffgf8WEZ+f5bm2A9sBent71x44cKCmzJOTk1xxxRVVHfPN8fP8xZPnOPfiy22Lfgne/88X8dYlC2vKUYtasmeFszdfu+YGZ2+Um2666fGI6J/ZPu8nAknfAK4ts+uOCp9bZdouqT6S7gAuAF+Y7UEiYj+wH6C/vz9qHSIpFApVD6/cseehS4oAwLkX4atPd/EHv1lbjlrUkj0rnL352jU3OHuzzVsIIuLts+2T9CNJiyPijKTFwDNluo0Cy0q2lwLjJY+xDXgncHOk+XjSQJ7Hb2adLO300cPAtuT+NuDLZfo8CqyUtELSImBrchySNlA8p3BLRJxNmaVhPI/fzNI6dHSMG/c8lMn1wtIWgj3AOySdpDgraA+ApCWSBgEi4gKwAxgCTgD3R8Tx5Pg/A14DPCDp25L2pczTEJ7Hb2ZpTK8XNjYxlcn1wlLNGoqI54Cby7SPA5tKtgeBwTL9/kma52+W6dlBnjVkZrWYa72wLLyPeK2hCnmtGrPa5X36ddbPM7oQWEXy/kK22vlSocXziWNl3vSzcp7Raw3ZvLI+vmnZ5mXUs3+e0YXA5uUXsqWR9WGRZsj6NTE8NGTz8gvZ0sj6sEizZPk8oz8R2Lz8PQpLI+vDIuZCYBXwC9nSyPqwiHloyCrg71FYWlkeFjEXAquQX8hmnctDQ2ZmOedCYGaWcy4EZmY550JgZpZzLgRmZjnnQmBmlnOpCoGkqyU9IOlk8vOqWfptkDQsaUTSrjL7f09SSLomTR4zM6te2k8Eu4AHI2Il8GCyfQlJXcBdwEZgNXCrpNUl+5dRvLrZ0ymzmJlZDdIWgs3Avcn9e4EtZfqsA0Yi4lREnAMOJMdN+1PgY0AmL1xvZtbp0haC3og4A5D8fH2ZPn3A6ZLt0aQNSbcAYxHxnZQ5zMysRvMuMSHpG8C1ZXbdUeFzqExbSLoseYxfr+hBpO3AdoDe3l4KhUKFT3+pycnJmo9tNWdvjXbN3q65wdmbLiJqvgHDwOLk/mJguEyftwBDJdu7k9ubgGeAHyS3CxTPE1w73/OuXbs2anXkyJGaj201Z2+Nds3errkjnL1RgMeizHtq2qGhw8C25P424Mtl+jwKrJS0QtIiYCtwOCKORcTrI2J5RCynOGR0fUT8MGUmMzOrQtpCsAd4h6STFGf+7AGQtETSIEBEXAB2AEPACeD+iDie8nnNzKxOUi1DHRHPATeXaR8HNpVsDwKD8zzW8jRZzMysNr4eQZMcOjrmC7uYWSa5EDTBoaNj7D54jKnzFwEYm5hi98FjAC4GZtZyXmuoCfYODb9UBKZNnb/I3qHhFiUyM3uZC0ETjE9MVdVuZtZMLgRNsKSnu6p2M7NmciFogp3rV9G9sOuStu6FXexcv6pFiczMXuaTxU0wfULYs4bMauNZd43lQtAkW9b0+RfXrAaeddd4Hhoys0zzrLvGcyEws0zzrLvGcyEws0zzrLvGcyEws0zzrLvG88liM8s0z7prPBcCM8s8z7prLA8NmZnlnAuBmVnOpSoEkq6W9ICkk8nPq2bpt0HSsKQRSbtm7PtQsu+4pE+myWNmZtVL+4lgF/BgRKwEHky2LyGpC7gL2AisBm6VtDrZdxOwGfgXEfHPgE+lzGNmZlVKWwg2A/cm9+8FtpTpsw4YiYhTEXEOOJAcB/DbwJ6IeAEgIp5JmcfMzKqUthD0RsQZgOTn68v06QNOl2yPJm0AbwT+taRHJP0PSb+aMo+ZmVVp3umjkr4BXFtm1x0VPofKtEXJ818F3AD8KnC/pDdERMw8QNJ2YDtAb28vhUKhwqe/1OTkZM3Htpqzt0a7Zm/X3ODsTRcRNd+AYWBxcn8xMFymz1uAoZLt3cDu5P7XgIGSfd8DXjff865duzZqdeTIkZqPbTVnb412zd6uuSOcvVGAx6LMe2raoaHDwLbk/jbgy2X6PAqslLRC0iJga3IcwCHgbQCS3ggsAn6cMpOZmVUhbSHYA7xD0kngHck2kpZIGgSIiAvADmAIOAHcHxHHk+PvAd4g6UmKJ5G3JVXLzMyaJNUSExHxHHBzmfZxYFPJ9iAwWKbfOeA9aTKYmVk6/maxmVnOuRCYmeWcC4GZWc65EJiZ5ZwLgZlZzrkQmJnlnAuBmVnOuRCYmeWcr1lsZtYGDh0dY+/QMOMTUyzp6Wbn+lV1u46zC4GZWcYdOjrG7oPHmDp/EYCxiSl2HzwGUJdi4KEhM7OM2zs0/FIRmDZ1/iJ7h4br8vguBGZmGTc+MVVVe7VcCMzMMm5JT3dV7dVyITAzy7id61fRvbDrkrbuhV3sXL+qLo/vk8VmZhk3fULYs4bMzHJsy5q+ur3xz5RqaEjS1ZIekHQy+XnVLP02SBqWNCJpV0n7dZIelvRtSY9JWpcmj5mZVS/tOYJdwIMRsRJ4MNm+hKQu4C5gI7AauFXS6mT3J4E/jojrgD9Mts3MrInSFoLNwL3J/XuBLWX6rANGIuJUcmnKA8lxAAFcmdz/R8B4yjxmZlYlpblWvKSJiOgp2f5JRFw1o8+7gQ0R8YFk+73AmyNih6RfoXhRe1EsSm+NiP83y3NtB7YD9Pb2rj1w4EBNmScnJ7niiitqOrbVnL012jV7u+YGZ2+Um2666fGI6J/ZPu/JYknfAK4ts+uOCp9bZdqmq89vAx+NiC9J+vfA54C3l3uQiNgP7Afo7++PgYGBCp/+UoVCgVqPbTVnb412zd6uucHZm23eQhARZd+YAST9SNLiiDgjaTHwTJluo8Cyku2lvDwEtA34cHL/r4G7K0ptZmZ1k/YcwWGKb+YkP79cps+jwEpJKyQtArYmx0GxIPyb5P7bgJMp85iZWZXSfo9gD3C/pN8Cngb+HYCkJcDdEbEpIi5I2kHxXEAXcE9EHE+O/0/ApyUtAH5Bcg7AzC7VyCWIzVIVgoh4Dri5TPs4sKlkexAYLNPv74C1aTKYdbpGL0Fs5rWGzDKu0UsQm7kQmGVco5cgNnMhMMu4Ri9BbOZCYJZxjV6C2Myrj5plXKOXIDZzITBrA41cgtjMQ0NmZjnnQmBmlnMuBGZmOedCYGaWcy4EZmY5l+rCNK0i6Vmg7AVsKnAN8OM6xmkmZ2+Nds3errnB2RvllyPidTMb27IQpCHpsXJX6GkHzt4a7Zq9XXODszebh4bMzHLOhcDMLOfyWAj2tzpACs7eGu2avV1zg7M3Ve7OEZiZ2aXy+InAzMxKdFQhkLRM0hFJJyQdl/ThpP1qSQ9IOpn8vKrkmN2SRiQNS1rfuvQv5emSdFTSV5LttsguqUfSFyV9N/n/f0sbZf9o8vvypKS/kvTqrGaXdI+kZyQ9WdJWdVZJayUdS/Z9RpJalH1v8jvzhKS/kdSTtezlcpfs+z1JIemarOWuSkR0zA1YDFyf3H8N8PfAauCTwK6kfRfwieT+auA7wKuAFcD3gK4W/xt+B/ivwFeS7bbIDtwLfCC5vwjoaYfsQB/wfaA72b4feH9WswO/BlwPPFnSVnVW4P8CbwEE/C2wsUXZfx1YkNz/RBazl8udtC8Dhih+p+marOWu5tZRnwgi4kxEfCu5/zPgBMUX+maKb1QkP7ck9zcDByLihYj4PjACrGtu6pdJWgr8BnB3SXPms0u6kuKL5XMAEXEuIiZog+yJBUC3pAXAZcA4Gc0eEf8T+IcZzVVllbQYuDIi/k8U36HuKzmmqdkj4usRcSHZfBhYmrXss/yfA/wp8DGg9ERrZnJXo6MKQSlJy4E1wCNAb0ScgWKxAF6fdOsDTpccNpq0tcqdFH+xXixpa4fsbwCeBf48Gda6W9LltEH2iBgDPgU8DZwBno+Ir9MG2UtUm7UvuT+zvdX+I8W/lCHj2SXdAoxFxHdm7Mp07tl0ZCGQdAXwJeAjEfHTubqWaWvJNCpJ7wSeiYjHKz2kTFurpoAtoPjR+bMRsQb4OcUhitlkJnsynr6Z4sf4JcDlkt4z1yFl2rI69W62rJn7N0i6A7gAfGG6qUy3TGSXdBlwB/CH5XaXactE7rl0XCGQtJBiEfhCRBxMmn+UfDQj+flM0j5KcZxv2lKKwwKtcCNwi6QfAAeAt0n6PO2RfRQYjYhHku0vUiwM7ZD97cD3I+LZiDgPHATeSntkn1Zt1lFeHoIpbW8JSduAdwL/IRk2gWxn/8cU/3D4TvJ6XQp8S9K1ZDv3rDqqECRn4T8HnIiIPynZdRjYltzfBny5pH2rpFdJWgGspHhCp+kiYndELI2I5cBW4KGIeA/tkf2HwGlJ01dTvxl4ijbITnFI6AZJlyW/PzdTPLfUDtmnVZU1GT76maQbkn/z+0qOaSpJG4DfB26JiLMluzKbPSKORcTrI2J58nodpThJ5YdZzj2nVp+trucN+FcUP249AXw7uW0CXgs8CJxMfl5dcswdFM/sD5ORs/jAAC/PGmqL7MB1wGPJ//0h4Ko2yv7HwHeBJ4G/pDjjI5PZgb+ieC7jPMU3oN+qJSvQn/x7vwf8GcmXS1uQfYTimPr063Vf1rKXyz1j/w9IZg1lKXc1N3+z2Mws5zpqaMjMzKrnQmBmlnMuBGZmOedCYGaWcy4EZmY550JgZpZzLgRmZjnnQmBmlnP/Hz7YAnz5I7lzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(nu,yn-yobs,'o')\n",
    "plt.plot([nu[0],nu[-1]],[0,0],'k-')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAV+klEQVR4nO3df4zcd33n8ee7jnO4LLHpOd3kbLeGO7cSjSWIt0l63J1mKa0SJyL3R3QKCqTJCa0SAYKeo4tbdHCcdLoUml6JgmJZwDU5OFaUhCpyjArNZUvyhwO7qZONMVwN9V3t+BygwmGJBdrL+/6Yr9tlmd39zsx3dnY/PB/SyDPz/cz3+/p+9+vXzHznV2QmkqQy/dywA0iSBseSl6SCWfKSVDBLXpIKZslLUsEuGtaCt27dmjt37hzoMn74wx/y6le/eqDL6IW5umOu7pirO+st18zMzHcz89LaM8rMoZz27NmTg/bEE08MfBm9MFd3zNUdc3VnveUCprOLrvVwjSQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSpY7ZKPiA0R8VcRcajDtIiI+yLiREQ8FxFXNhtTktSLbh7Jvw84vsS064Bd1WkCeKDPXJKkBtQq+YjYDlwPfGKJITcCD1Xv1T8CbImIyxvKKEnqUWSNHw2JiM8D/wV4DXBXZt6waPoh4J7MfKq6/Dhwd2ZOLxo3QfuRPqOjo3smJycbWYmlzM3NMTIyMtBl9GK95Zo9fW4IaWD3ts3AcLZXnXUe3QRnzze73Avr3I/1tn8N23rLNT4+PpOZY3Xns+J310TEDcCLmTkTEa2lhnW47qfuPTLzIHAQYGxsLFutpWbXjKmpKQa9jF6st1y37X9s9cMAJ29pAcPZXnXWed/uee6dbfbrny6scz/W2/41bKXnqnO45s3A2yLiJDAJvCUiPr1ozClgx4LL24EX+k4nSerLiiWfmb+XmdszcydwM/A/M/Mdi4Y9CtxavcvmGuBcZp5pPq4kqRs9P9eMiDsAMvMAcBjYC5wAXgZubySdJKkvXZV8Zk4BU9X5AwuuT+DdTQaTJPXPT7xKUsEseUkqmCUvSQWz5CWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgq2YslHxKsi4qsR8WxEHIuID3cY04qIcxFxtDp9cDBxJUndqPPLUD8C3pKZcxGxEXgqIr6YmUcWjXsyM29oPqIkqVcrlnz1035z1cWN1SkHGUqS1Ixax+QjYkNEHAVeBL6cmU93GPYb1SGdL0bErzWaUpLUk2g/UK85OGIL8AXgvZn5/ILrLwFeqQ7p7AU+lpm7Otx+ApgAGB0d3TM5Odlv/mXNzc0xMjIy0GX0Yr3lmj19bghpYPe2zcBwtleddR7dBGfPN7vcC+vcj/W2fw3bess1Pj4+k5ljdefTVckDRMSHgB9m5h8uM+YkMJaZ311qzNjYWE5PT3e17G5NTU3RarUGuoxerLdcO/c/tvphgJP3XA8MZ3vVWed9u+e5d7bOy1r1XVjnfqy3/WvY1luuiOiq5Ou8u+bS6hE8EbEJeCvwjUVjLouIqM5fVc33e3VDSJIGo87DkMuBByNiA+3y/lxmHoqIOwAy8wBwE3BnRMwD54Gbs9unCJKkxtV5d81zwJs6XH9gwfn7gfubjSZJ6pefeJWkglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SC1fmN11dFxFcj4tmIOBYRH+4wJiLivog4ERHPRcSVg4krSepGnd94/RHwlsyci4iNwFMR8cXMPLJgzHXArup0NfBA9a8kaYhWfCSfbXPVxY3VafGPdN8IPFSNPQJsiYjLm40qSepWZC7u6w6DIjYAM8A/Az6emXcvmn4IuCczn6ouPw7cnZnTi8ZNABMAo6OjeyYnJxtZiaXMzc0xMjIy0GX0Yr3lmj19bghpYPe2zcBwtleddR7dBGfPN7vcC+vcj/W2fw3bess1Pj4+k5ljdedT53ANmfn/gDdGxBbgCxFxRWY+v2BIdLpZh/kcBA4CjI2NZavVqpuzJ1NTUwx6Gb1Yb7lu2//Y6ocBTt7SAoazveqs877d89w7W+u/UG0X1rkf623/GrbSc3X17prM/D4wBVy7aNIpYMeCy9uBF/pKJknqW51311xaPYInIjYBbwW+sWjYo8Ct1btsrgHOZeaZxtNKkrpS57nm5cCD1XH5nwM+l5mHIuIOgMw8ABwG9gIngJeB2weUV5LUhRVLPjOfA97U4foDC84n8O5mo0mS+uUnXiWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalgdX7jdUdEPBERxyPiWES8r8OYVkSci4ij1emDg4krSepGnd94nQf2ZeYzEfEaYCYivpyZX1807snMvKH5iJKkXq34SD4zz2TmM9X5HwDHgW2DDiZJ6l+0f4O75uCIncBXgCsy86UF17eAh4FTwAvAXZl5rMPtJ4AJgNHR0T2Tk5N9RF/Z3NwcIyMjA11GL9ZbrtnT54aQBnZv2wwMZ3vVWefRTXD2fLPLvbDO/Vhv+9ewrbdc4+PjM5k5Vnc+tUs+IkaAvwT+c2Y+smjaJcArmTkXEXuBj2XmruXmNzY2ltPT03Vz9mRqaopWqzXQZfRiveXauf+x1Q8DnLznemA426vOOu/bPc+9s3WOeNZ3YZ37sd72r2Fbb7kioquSr/XumojYSPuR+mcWFzxAZr6UmXPV+cPAxojYWjeEJGkw6ry7JoBPAscz84+WGHNZNY6IuKqa7/eaDCpJ6l6d55pvBt4JzEbE0eq63wd+CSAzDwA3AXdGxDxwHrg5uznYL0kaiBVLPjOfAmKFMfcD9zcVSpLUDD/xKkkFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SCWfKSVDBLXpIKZslLUsEseUkqmCUvSQWr8xuvOyLiiYg4HhHHIuJ9HcZERNwXESci4rmIuHIwcSVJ3ajzG6/zwL7MfCYiXgPMRMSXM/PrC8ZcB+yqTlcDD1T/SpKGaMVH8pl5JjOfqc7/ADgObFs07EbgoWw7AmyJiMsbTytJ6kpkZv3BETuBrwBXZOZLC64/BNxT/eg3EfE4cHdmTi+6/QQwATA6OrpncnKyp9Czp8/VGje6Cc6e72kRS9q9bXPf85ibm2NkZKSr29Rd534MYns14Wcp17D2r9WwVK7V2LeXM4z9q87feantNT4+PpOZY3WXVedwDQARMQI8DLx/YcFfmNzhJj9175GZB4GDAGNjY9lqteou/ifctv+xWuP27Z7n3tnaq1jLyVtafc9jamqKbte97jr3YxDbqwk/S7mGtX+thqVyrca+vZxh7F91/s5N/R1rvbsmIjbSLvjPZOYjHYacAnYsuLwdeKHvdJKkvtR5d00AnwSOZ+YfLTHsUeDW6l021wDnMvNMgzklST2o8xzlzcA7gdmIOFpd9/vALwFk5gHgMLAXOAG8DNzefFRJUrdWLPnqxdROx9wXjkng3U2FkiQ1w0+8SlLBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SCWfKSVDBLXpIKZslLUsHq/MbrpyLixYh4fonprYg4FxFHq9MHm48pSepFnd94/RPgfuChZcY8mZk3NJJIktSYFR/JZ+ZXgL9bhSySpIZF+ze4VxgUsRM4lJlXdJjWAh4GTgEvAHdl5rEl5jMBTACMjo7umZyc7Cn07OlztcaNboKz53taxJJ2b9vc9zzm5uYYGRnp6jZ117kfg9heTfhZyjWs/Ws1LJVrNfbt5Qxj/6rzd15qe42Pj89k5ljdZTVR8pcAr2TmXETsBT6WmbtWmufY2FhOT0/XzfkTdu5/rNa4fbvnuXe2zhGp+k7ec33f85iamqLVanV1m7rr3I9BbK8m/CzlGtb+tRqWyrUa+/ZyhrF/1fk7L7W9IqKrku/73TWZ+VJmzlXnDwMbI2Jrv/OVJPWv75KPiMsiIqrzV1Xz/F6/85Uk9W/F5ygR8VmgBWyNiFPAh4CNAJl5ALgJuDMi5oHzwM1Z5xiQJGngViz5zHz7CtPvp/0WS0nSGuMnXiWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalgK5Z8RHwqIl6MiOeXmB4RcV9EnIiI5yLiyuZjSpJ6UeeR/J8A1y4z/TpgV3WaAB7oP5YkqQkrlnxmfgX4u2WG3Ag8lG1HgC0RcXlTASVJvYvMXHlQxE7gUGZe0WHaIeCezHyquvw4cHdmTncYO0H70T6jo6N7Jicnewo9e/pcrXGjm+Ds+Z4WsaTd2zb3PY+5uTlGRka6uk3dde7HILZXE8zVHXN1Zxi56vTIUj0xPj4+k5ljdZd1UXfROooO13W858jMg8BBgLGxsWy1Wj0t8Lb9j9Uat2/3PPfONrGK/+DkLa2+5zE1NUW36153nfsxiO3VBHN1x1zdGUauOj3SS0900sS7a04BOxZc3g680MB8JUl9aqLkHwVurd5lcw1wLjPPNDBfSVKfVnyOEhGfBVrA1og4BXwI2AiQmQeAw8Be4ATwMnD7oMJKkrqzYsln5ttXmJ7AuxtLJElqjJ94laSCWfKSVDBLXpIKZslLUsEseUkqmCUvSQWz5CWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpILVKvmIuDYivhkRJyJif4fprYg4FxFHq9MHm48qSepWnd943QB8HPgt4BTwtYh4NDO/vmjok5l5wwAySpJ6VOeR/FXAicz8dmb+GJgEbhxsLElSE6L9O9zLDIi4Cbg2M99VXX4ncHVmvmfBmBbwMO1H+i8Ad2XmsQ7zmgAmAEZHR/dMTk72FHr29Lla40Y3wdnzPS1iSbu3be57HnNzc4yMjHR1m7rr3I9BbK8mmKs75urOMHLV6ZGlemJ8fHwmM8fqLmvFwzVAdLhu8T3DM8AvZ+ZcROwF/gzY9VM3yjwIHAQYGxvLVqtVN+dPuG3/Y7XG7ds9z72zdVaxvpO3tPqex9TUFN2ue9117scgtlcTzNUdc3VnGLnq9EgvPdFJncM1p4AdCy5vp/1o/e9l5kuZOVedPwxsjIitfaeTJPWlTsl/DdgVEa+LiIuBm4FHFw6IiMsiIqrzV1Xz/V7TYSVJ3VnxOUpmzkfEe4A/BzYAn8rMYxFxRzX9AHATcGdEzAPngZtzpYP9kqSBq3UgqjoEc3jRdQcWnL8fuL/ZaJKkfvmJV0kqmCUvSQWz5CWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSpYrZKPiGsj4psRcSIi9neYHhFxXzX9uYi4svmokqRurVjyEbEB+DhwHfAG4O0R8YZFw64DdlWnCeCBhnNKknpQ55H8VcCJzPx2Zv4YmARuXDTmRuChbDsCbImIyxvOKknqUmTm8gMibgKuzcx3VZffCVydme9ZMOYQcE9mPlVdfhy4OzOnF81rgvYjfYBfBb7Z1IosYSvw3QEvoxfm6o65umOu7qy3XL+cmZfWnclFNcZEh+sW3zPUGUNmHgQO1lhmIyJiOjPHVmt5dZmrO+bqjrm6U3quOodrTgE7FlzeDrzQwxhJ0iqrU/JfA3ZFxOsi4mLgZuDRRWMeBW6t3mVzDXAuM880nFWS1KUVD9dk5nxEvAf4c2AD8KnMPBYRd1TTDwCHgb3ACeBl4PbBRe7Kqh0a6pK5umOu7pirO0XnWvGFV0nS+uUnXiWpYJa8JBVs3Zd8RPxCRHw5Iv66+ve1S4zr+NUMEfHGiDgSEUcjYjoirloLuapp762mHYuIj6yVXNX0uyIiI2LrWsgVER+NiG9UX6vxhYjY0meenr/KY6XbDiNXROyIiCci4ni1P71vLeRaMH1DRPxV9ZmbNZErIrZExOer/ep4RPzGGsn1u9Xf8PmI+GxEvGrZhWXmuj4BHwH2V+f3A3/QYcwG4FvA64GLgWeBN1TTvgRcV53fC0ytkVzjwF8A/6i6/ItrIVc1fQftF+L/N7B1LeQCfhu4qDr/B51u30WWZdd/wb7yRdqfEbkGeLrubYeU63Lgyur8a4D/tRZyLZj+74D/ARxqIlMTuYAHgXdV5y8Gtgw7F7AN+BtgU3X5c8Btyy1v3T+Sp/2VCg9W5x8E/nWHMct9NUMCl1TnN9Pc+/v7zXUn7U8R/wggM19cI7kA/ivw7+nwgbdh5crML2XmfDXuCO3PavSqn6/yqHPbVc+VmWcy8xmAzPwBcJx2YQw1F0BEbAeuBz7RUJ6+c0XEJcC/Aj4JkJk/zszvDztXNe0iYFNEXAT8PCt0VgklP5rVe/Krf3+xw5htwN8uuHyKf9jB3w98NCL+FvhD4PfWSK5fAf5lRDwdEX8ZEb++FnJFxNuA05n5bEN5Gsm1yL+l/SioV3WWs9SYuhlXO9ffi4idwJuAp9dIrj+m/aDhlYbyNJHr9cB3gP9WHUb6RES8eti5MvM07Z76P8AZ2p9J+tJyC6vztQZDFxF/AVzWYdIH6s6iw3UXHoXeCfxuZj4cEf+G9j33W9dArouA19J+qvbrwOci4vVZPUcbRq6I+PlqHr9dcz6rkmvRMj4AzAOf6S5dd8tZZkytr/joUd9fMRIRI8DDwPsz86Vh54qIG4AXM3MmIloN5ek7F+3/f1cC783MpyPiY7QPI/6HYeaqXqu6EXgd8H3gTyPiHZn56aUWti5KPjOXLN2IOHvh6Wj1dKbTYY3lvnbhd4ALL0L9KV08ZRxwrlPAI1WpfzUiXqH9hUXfGWKuf0p753o2Ii5c/0xEXJWZ/3eIuS7M43eAG4DfrHNnuIx+vsrj4hq3HUYuImIj7YL/TGY+0lCmfnPdBLwtIvYCrwIuiYhPZ+Y7hpwrgVOZeeHZzudpl3wT+sn1VuBvMvM7ABHxCPDPgSVLvpEXOIZ5Aj7KT75g95EOYy4Cvk27oC680PFr1bTjQKs6/5vAzBrJdQfwn6rzv0L7qVsMO9eicSdp7oXXfrfXtcDXgUsbyLLi+tM+hrzwhbGvdrPthpArgIeAP24iS1O5Fo1p0ewLr33lAp4EfrU6/x+Bjw47F3A1cIz2sfig/frVe5ddXtN/8NU+Af8YeBz46+rfX6iu/yfA4QXj9tJ+R8G3gA8suP5fADPVhn4a2LNGcl1M+975eeAZ4C1rIdeieZ2kuZLvd3udoH1HeLQ6Hegzz08th/Yd7x3V+aD9YzrfAmaBsW623WrnqvbzBJ5bsI32DjvXonm0aLDkG/g7vhGYrrbZnwGvXSO5Pgx8g3Y3/Heqd+AtdfJrDSSpYCW8u0aStARLXpIKZslLUsEseUkqmCUvSQWz5CWpYJa8JBXs/wNBnfpOsvYdUAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(yn-yobs)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
